! This module is useful for verification tests and development of
! turbulence models.

MODULE TURBULENCE

USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_POINTERS
USE MESH_VARIABLES
USE COMP_FUNCTIONS

IMPLICIT NONE

CHARACTER(255), PARAMETER :: turbid='$Id: turb.f90 2900 2008-12-18 18:26:05Z drjfloyd $'
CHARACTER(255), PARAMETER :: turbrev='$Revision: 2900 $'
CHARACTER(255), PARAMETER :: turbdate='$Date: 2008-12-18 13:26:05 -0500 (Thu, 18 Dec 2008) $'

PRIVATE
PUBLIC :: ANALYTICAL_SOLUTION, sandia_dat, INIT_TURB_ARRAYS, spectral_output, VARDEN_DYNSMAG, &
          GET_REV_turb, TURBULENT_KINETIC_ENERGY, WERNER_WENGLE_WALL_MODEL, COMPRESSION_WAVE, &
          DUMP_RAW_DATA, SURFACE_HEAT_FLUX_MODEL, SUBGRID_SCALAR_VARIANCE, WAVELET_ERROR_ESTIMATOR
 
CONTAINS


SUBROUTINE SANDIA_DAT(NM)
IMPLICIT NONE

! This routine reads the file 'iso_ini.dat', which is generated by turb_init.
! This exe generates a random velocity field with a spectrum that matches the
! Comte-Bellot/Corrsin 1971 experimental data.

REAL :: XXX
INTEGER :: I,J,K,II,JJ,KK,FILE_NUM
INTEGER, INTENT(IN) :: NM
LOGICAL :: CBC=.FALSE.

CALL POINT_TO_MESH(NM)

IF (XF<1._EB) CBC=.TRUE. ! XF=9*2*pi/100 meters for CBC

IF (IBAR/=32 .AND. IBAR/=64) THEN
   WRITE(LU_ERR,'(A)') 'Error 1 in SANDIA_DAT!'
   STOP
ENDIF  

FILE_NUM = GET_FILE_NUMBER()
IF (CBC) THEN
   IF (IBAR==32) OPEN (UNIT=FILE_NUM,FILE='iso_ini.32.cbc.rjm.dat',FORM='formatted',STATUS='old')
   IF (IBAR==64) OPEN (UNIT=FILE_NUM,FILE='iso_ini.64.cbc.rjm.dat',FORM='formatted',STATUS='old')
ELSE
   IF (IBAR==32) OPEN (UNIT=FILE_NUM,FILE='iso_ini.32.kang.rjm.dat',FORM='formatted',STATUS='old')
   IF (IBAR==64) OPEN (UNIT=FILE_NUM,FILE='iso_ini.64.kang.rjm.dat',FORM='formatted',STATUS='old') 
ENDIF
     
READ (FILE_NUM,*) II, JJ, KK	! reads number of points in each direction

IF (II/=IBAR .OR. JJ/=JBAR .OR. KK/=KBAR) THEN
   WRITE(LU_ERR,'(A)') 'Error 2 in SANDIA_DAT!'
   STOP
ENDIF

READ (FILE_NUM,*) XXX, XXX, XXX	! reads lower physical dimension limit
READ (FILE_NUM,*) XXX, XXX, XXX	! reads upper physical dimension limit

DO K = 1,KBAR
   DO J = 1,JBAR
      DO I = 1,IBAR

         READ (FILE_NUM,*) U(I,J,K), V(I,J,K), W(I,J,K), XXX, XXX
         
         IF (CBC) THEN
            U(I,J,K) = 0.01*U(I,J,K)	! scale from cm/s
            V(I,J,K) = 0.01*V(I,J,K)
            W(I,J,K) = 0.01*W(I,J,K)
         ENDIF
         
      ENDDO
   ENDDO
ENDDO

CLOSE (UNIT=FILE_NUM)

! subtract mean
U(1:IBAR,1:JBAR,1:KBAR) = U(1:IBAR,1:JBAR,1:KBAR) - SUM(U(1:IBAR,1:JBAR,1:KBAR))/REAL(IBAR*JBAR*KBAR,EB)
V(1:IBAR,1:JBAR,1:KBAR) = V(1:IBAR,1:JBAR,1:KBAR) - SUM(V(1:IBAR,1:JBAR,1:KBAR))/REAL(IBAR*JBAR*KBAR,EB)
W(1:IBAR,1:JBAR,1:KBAR) = W(1:IBAR,1:JBAR,1:KBAR) - SUM(W(1:IBAR,1:JBAR,1:KBAR))/REAL(IBAR*JBAR*KBAR,EB)

! apply periodic b.c.
U(0,:,:) = U(IBAR,:,:)
V(:,0,:) = V(:,JBAR,:)
W(:,:,0) = W(:,:,KBAR)

END SUBROUTINE SANDIA_DAT


SUBROUTINE INIT_TURB_ARRAYS(NM)
USE MEMORY_FUNCTIONS, ONLY: ChkMemErr
IMPLICIT NONE
INTEGER, INTENT(IN) :: NM
INTEGER :: IZERO
TYPE (MESH_TYPE), POINTER :: M

CALL POINT_TO_MESH(NM)
M => MESHES(NM)

IF (PERIODIC_TEST==2 .OR. DYNSMAG) THEN
   ! real work arrays
   ALLOCATE(M%TURB_WORK1(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK1',IZERO)
   M%TURB_WORK1 = 0._EB
   ALLOCATE(M%TURB_WORK2(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK2',IZERO)
   M%TURB_WORK2 = 0._EB
   ALLOCATE(M%TURB_WORK3(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK3',IZERO)
   M%TURB_WORK3 = 0._EB
   ALLOCATE(M%TURB_WORK4(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK4',IZERO)
   M%TURB_WORK4 = 0._EB
ENDIF

IF (PERIODIC_TEST==2) THEN
   ! complex work arrays
   ALLOCATE(M%TURB_WORK_DPC1(IBAR,JBAR,KBAR),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK_DPC1',IZERO)
   M%TURB_WORK_DPC1 = 0._EB
   ALLOCATE(M%TURB_WORK_DPC2(IBAR,JBAR,KBAR),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK_DPC2',IZERO)
   M%TURB_WORK_DPC2 = 0._EB
   ALLOCATE(M%TURB_WORK_DPC3(IBAR,JBAR,KBAR),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK_DPC3',IZERO)
   M%TURB_WORK_DPC3 = 0._EB
   ALLOCATE(M%TURB_WORK_DPC4(IBAR,JBAR,KBAR),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK_DPC4',IZERO)
   M%TURB_WORK_DPC4 = 0._EB
ENDIF

IF (DYNSMAG) THEN
   ! real work arrays
   ALLOCATE(M%TURB_WORK5(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK5',IZERO)
   M%TURB_WORK5 = 0._EB
   ALLOCATE(M%TURB_WORK6(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK6',IZERO)
   M%TURB_WORK6 = 0._EB
   ALLOCATE(M%TURB_WORK7(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK7',IZERO)
   M%TURB_WORK7 = 0._EB
   ALLOCATE(M%TURB_WORK8(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK8',IZERO)
   M%TURB_WORK8 = 0._EB
   
   ALLOCATE(M%TURB_WORK9(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK9',IZERO)
   M%TURB_WORK9 = 0._EB
   ALLOCATE(M%TURB_WORK10(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK10',IZERO)
   M%TURB_WORK10 = 0._EB
ENDIF
   
! 1D working arrays
IF (DYNSMAG .OR. CHECK_KINETIC_ENERGY) THEN
   ALLOCATE(M%TURB_WORK11(MAX(IBAR,JBAR,KBAR)),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK11',IZERO)
   M%TURB_WORK11 = 0._EB
   ALLOCATE(M%TURB_WORK12(MAX(IBAR,JBAR,KBAR)),STAT=IZERO)
   CALL ChkMemErr('INIT_TURB_ARRAYS','TURB_WORK12',IZERO)
   M%TURB_WORK12 = 0._EB
ENDIF

END SUBROUTINE INIT_TURB_ARRAYS


SUBROUTINE SPECTRAL_OUTPUT(NM)
IMPLICIT NONE
INTEGER, INTENT(IN) :: NM
INTEGER :: nn(3),N_LO(3),N_HI(3),ARRAY_LO(3),ARRAY_HI(3)
REAL(EB),     POINTER, DIMENSION(:,:,:) :: UU,VV,WW,HH
COMPLEX(DPC), POINTER, DIMENSION(:,:,:) :: UUHT,VVHT,WWHT,KKHT

call POINT_TO_MESH(NM)

nn = IBAR

N_LO = 1
N_HI = (/IBAR,JBAR,KBAR/)

ARRAY_LO = 0
ARRAY_HI = (/IBP1,JBP1,KBP1/)

UU => TURB_WORK1
VV => TURB_WORK2
WW => TURB_WORK3
HH => TURB_WORK4
UUHT => TURB_WORK_DPC1
VVHT => TURB_WORK_DPC2
WWHT => TURB_WORK_DPC3
KKHT => TURB_WORK_DPC4

UU = U
VV = V
WW = W
HH = H

! take fourier transform of velocities in 3d...
call fft3d_f90(UU,UUHT,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
call fft3d_f90(VV,VVHT,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
call fft3d_f90(WW,WWHT,N_LO,N_HI,ARRAY_LO,ARRAY_HI)

! calc the spectral kinetic energy
call complex_tke_f90(KKHT, UUHT, VVHT, WWHT, nn(1))

! total up the spectral energy for each mode and integrate over
! the resolved modes...
call spectrum_f90(KKHT, nn(1), XF-XS, nint(100._EB*SPEC_CLOCK))
      
IF (TURB_INIT) call sandia_out(UU,VV,WW,HH,nn(1))

IF (SPEC_CLOCK>=0.00_EB .AND. SPEC_CLOCK<0.28_EB) THEN
   SPEC_CLOCK=0.28_EB
   RETURN
ENDIF
IF (SPEC_CLOCK>=0.28_EB .AND. SPEC_CLOCK<0.66_EB) THEN
   SPEC_CLOCK=0.66_EB
   RETURN
ENDIF
IF (SPEC_CLOCK>=0.66_EB) SPEC_CLOCK=1.E9_EB

END SUBROUTINE SPECTRAL_OUTPUT


SUBROUTINE sandia_out(u,v,w,p,n)
IMPLICIT NONE

! Variable declarations

INTEGER, INTENT(IN) :: n
REAL(EB), INTENT(IN) :: u(n,n,n)
REAL(EB), INTENT(IN) :: v(n,n,n)
REAL(EB), INTENT(IN) :: w(n,n,n)
REAL(EB), INTENT(IN) :: p(n,n,n)
REAL(EB) :: uu,vv,ww,pp,ke

INTEGER :: i,j,k,file_num

! This subroutine writes out the velocity, relative pressure, and
! turbulent kinetic energy to the file 'ini_salsa.dat'.  This is
! then used as input to the 'turb_init' program.

file_num = GET_FILE_NUMBER()
OPEN (unit=file_num, file='ini_salsa.dat', form='formatted', status='unknown', position='rewind')

WRITE (file_num,997) n,n,n
WRITE (file_num,998) 0,0,0
WRITE (file_num,998) 2*pi,2*pi,2*pi

DO k = 1,n
   DO j = 1,n
      DO i = 1,n

         pp = 10._EB*p(i,j,k)  ! convert to dynes from pascals
         uu = 100._EB*u(i,j,k) ! convert to cm/s from m/s
         vv = 100._EB*v(i,j,k)
         ww = 100._EB*w(i,j,k)
         	
         ke = 0.5_EB*( uu**2 + vv**2 + ww**2 )

         WRITE (file_num,999) uu,vv,ww,pp,ke
         
      END DO
   END DO
END DO

CLOSE (unit=file_num)


997	FORMAT(3(i6,8x))
998	FORMAT(3(f12.6,2x))
999	FORMAT(5(f12.6,2x))


END SUBROUTINE sandia_out


SUBROUTINE complex_tke_f90(tkeht, upht, vpht, wpht, n)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
COMPLEX(DPC), INTENT(OUT) :: tkeht(n,n,n)
COMPLEX(DPC), INTENT(IN) :: upht(n,n,n),vpht(n,n,n),wpht(n,n,n)
INTEGER i,j,k

do k = 1,n
   do j = 1,n
      do i = 1,n
         tkeht(i,j,k) = 0.5_EB*( upht(i,j,k)*conjg(upht(i,j,k)) + &
                                 vpht(i,j,k)*conjg(vpht(i,j,k)) + &
                                 wpht(i,j,k)*conjg(wpht(i,j,k)) )
      end do
   end do
end do

END SUBROUTINE complex_tke_f90


SUBROUTINE fft3d_f90(v,vht,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
IMPLICIT NONE
INTEGER, INTENT(IN) :: N_LO(3),N_HI(3),ARRAY_LO(3),ARRAY_HI(3)
REAL(EB), INTENT(IN) :: v(ARRAY_LO(1):ARRAY_HI(1),ARRAY_LO(2):ARRAY_HI(2),ARRAY_LO(3):ARRAY_HI(3))
COMPLEX(DPC), INTENT(INOUT) :: vht(N_LO(1):N_HI(1),N_LO(2):N_HI(2),N_LO(3):N_HI(3))

! This routine performs an FFT on the real array v and places the
! result in the complex array vht.

INTEGER :: i,j,k,nn(3)
REAL(EB) :: z2(2*N_HI(1))
!REAL(EB) :: tke

nn = N_HI

! convert v to spectral space

!tke = 0._EB

do k = 1,nn(3)
   do j = 1,nn(2)
      do i = 1,nn(1)      
         vht(i,j,k) = cmplx(v(i,j,k),0.0_EB,kind=DPC)
         !tke = tke + 0.5_EB*v(i,j,k)**2
      end do
   end do
end do

!tke = tke/real(nn(1)*nn(2)*nn(3),EB)

!print*, ' fft3d internal check-'
!print*, ' tkeave (physical) =', tke

!do k=1,nn(3)
!   do j=1,nn(2)
!      do i=1,nn(1)
!         vht(i,j,k) = cmplx(float(i)+k*100,float(j)+.5+k*100)
!      enddo
!   enddo
!enddo

call fourt_f90(vht,nn,3,-1,0,z2)
!call four3(vht,-1)

do k = 1,nn(3)
   do j = 1,nn(2)
      do i = 1,nn(1)
         vht(i,j,k) = vht(i,j,k)/nn(1)/nn(2)/nn(3)
      end do
   end do
end do

!tke = 0._EB
!do k = 1,nn(3)
!   do j = 1,nn(2)
!      do i = 1,nn(1)
!         tke = tke + 0.5_EB*real(vht(i,j,k)*conjg(vht(i,j,k)),kind=EB)
!      end do
!   end do
!end do

!print*, ' tkeave (spectral) =', tke

END SUBROUTINE fft3d_f90


SUBROUTINE spectrum_f90(vht, n, Lm, iFVfilenum)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n, iFVfilenum 
COMPLEX(DPC), INTENT(IN) :: vht(n,n,n)
REAL(EB), INTENT(IN) :: Lm

! This routine is copied from SNL and is intended to compute the
! kinetic energy spectrum in wavenumber space.

INTEGER :: kmax, kx, ky, kz, k, ksum, file_num
INTEGER :: num(0:n)

REAL(EB) :: rk, temp, vsum, rkx, rky, rkz, etot
REAL(EB) :: vt(0:n)

CHARACTER*6 :: ext
CHARACTER*30 :: filename

! for dimensional wavenumbers
REAL(EB) :: wn(0:n)
REAL(EB) :: L, k0
      
L = Lm

k0 = 2*PI/L
kmax = n/2
	
wn(0) = 0._EB
do k = 1,n
   wn(k) = k0*k
end do

do k = 0,n
   vt(k) = 0._EB
   num(k) = 0
end do

etot = 0._EB

do kx = 1,n

   rkx = real(kx-1,kind=EB)
   if (rkx .gt. kmax) then
      rkx = n - rkx
   end if

   do ky = 1,n

      rky = real(ky-1,kind=EB)
      if (rky .gt. kmax) then
         rky = n - rky
      end if

      do kz = 1,n

         rkz = real(kz-1,kind=EB)
         if (rkz .gt. kmax) then
            rkz = n - rkz
         end if

         rk     = sqrt(rkx*rkx + rky*rky + rkz*rkz)
         k      = nint(rk)

         num(k) = num(k) + 1
         temp   = real(vht(kx,ky,kz)*conjg(vht(kx,ky,kz)),kind=EB)
         etot   = etot + sqrt(temp)
         vt(k)  = vt(k) + sqrt(temp)*(L/(2*PI))

      end do
   end do
end do

write(6,*) ' '
write(6,*) ' Spectrum Internal Check-'
write(6,*) ' Total Energy (-) in 3D field = ', etot
write(6,*) ' k(-), Num(k), k(1/cm), E(cm3/s2) '
ksum = 0
vsum = 0._EB
do k = 0,n
   write(6,*) k, num(k), wn(k), vt(k)
   ksum = ksum + num(k)
   vsum = vsum + vt(k)
end do

write(6,*) ' ksum: ', ksum
write(6,*) ' Total Energy (-) in spectrum: ', vsum
write(6,*) ' '

! write the spectral data to a file
! units are (1/cm) for wavenumber and (cm3/s2) for energy
! this matches the Comte-Bellot/Corrsin units

write(ext,1) iFVfilenum
filename = TRIM(CHID) // '_spec_' // TRIM(ext) // '.dat'

file_num = GET_FILE_NUMBER()
open (unit=file_num, file=filename, status='unknown', form='formatted')

do k = 0,n
   write (file_num,*) wn(k), vt(k)
end do

close (unit = file_num)

1	format (i3.3)

END SUBROUTINE spectrum_f90


SUBROUTINE fourt_f90(data3,nn,ndim,isign,iform,work)
IMPLICIT NONE
! Converted to F90 10/8/2008 by RJM
!
!     the cooley-tukey fast fourier transform in usasi basic fortran
!     transform(j1,j2,,,,) = sum(data(i1,i2,,,,)*w1**((i2-1)*(j2-1))
!                                 *w2**((i2-1)*(j2-1))*,,,),
!     where i1 and j1 run from 1 to nn(1) and w1=exp(isign*2*pi=
!     sqrt(-1)/nn(1)), etc.  there is no limit on the dimensionality
!     (number of subscripts) of the data array.  if an inverse
!     transform (isign=+1) is performed upon an array of transformed
!     (isign=-1) data, the original data will reappear.
!     multiplied by nn(1)*nn(2)*,,,  the array of input data must be
!     in complex format.  however, if all imaginary parts are zero (i.e.
!     the data are disguised real) running time is cut up to forty per-
!     cent.  (for fastest transform of real data, nn(1) should be even.)
!     the transform values are always complex and are returned in the
!     original array of data, replacing the input data.  the length
!     of each dimension of the data array may be any integer.  the
!     program runs faster on composite integers than on primes, and is
!     particularly fast on numbers rich in factors of two.
!
!     timing is in fact given by the following formula.  let ntot be the
!     total number of points (real or complex) in the data array, that
!     is, ntot=nn(1)*nn(2)*...  decompose ntot into its prime factors,
!     such as 2**k2 * 3**k3 * 5**k5 * ...  let sum2 be the sum of all
!     the factors of two in ntot, that is, sum2 = 2*k2.  let sumf be
!     the sum of all other factors of ntot, that is, sumf = 3*k3*5*k5*..
!     the time taken by a multidimensional transform on these ntot data
!     is t = t0 + ntot*(t1+t2*sum2+t3*sumf).  on the cdc 3300 (floating
!     point add time = six microseconds), t = 3000 + ntot*(600+40*sum2+
!     175*sumf) microseconds on complex data.
!
!     implementation of the definition by summation will run in a time
!     proportional to ntot*(nn(1)+nn(2)+...).  for highly composite ntot
!     the savings offered by this program can be dramatic.  a one-dimen-
!     sional array 4000 in length will be transformed in 4000*(600+
!     40*(2+2+2+2+2)+175*(5+5+5)) = 14.5 seconds versus about 4000*
!     4000*175 = 2800 seconds for the straightforward technique.
!
!     the fast fourier transform places three restrictions upon the
!     data.
!     1.  the number of input data and the number of transform values
!     must be the same.
!     2.  both the input data and the transform values must represent
!     equispaced points in their respective domains of time and
!     frequency.  calling these spacings deltat and deltaf, it must be
!     true that deltaf=2*pi/(nn(i)*deltat).  of course, deltat need not
!     be the same for every dimension.
!     3.  conceptually at least, the input data and the transform output
!     represent single cycles of periodic functions.
!
!     the calling sequence is--
!     call fourt(data,nn,ndim,isign,iform,work)
!
!     data is the array used to hold the real and imaginary parts
!     of the data on input and the transform values on output.  it
!     is a multidimensional floating point array, with the real and
!     imaginary parts of a datum stored immediately adjacent in storage
!     (such as fortran iv places them).  normal fortran ordering is
!     expected, the first subscript changing fastest.  the dimensions
!     are given in the integer array nn, of length ndim.  isign is -1
!     to indicate a forward transform (exponential sign is -) and +1
!     for an inverse transform (sign is +).  iform is +1 if the data are
!     complex, 0 if the data are real.  if it is 0, the imaginary
!     parts of the data must be set to zero.  as explained above, the
!     transform values are always complex and are stored in array data.
!     work is an array used for working storage.  it is floating point
!     real, one dimensional of length equal to twice the largest array
!     dimension nn(i) that is not a power of two.  if all nn(i) are
!     powers of two, it is not needed and may be replaced by zero in the
!     calling sequence.  thus, for a one-dimensional array, nn(1) odd,
!     work occupies as many storage locations as data.  if supplied,
!     work must not be the same array as data.  all subscripts of all
!     arrays begin at one.
!
!     example 1.  three-dimensional forward fourier transform of a
!     complex array dimensioned 32 by 25 by 13 in fortran iv.
!     dimension data(32,25,13),work(50),nn(3)
!     complex data
!     data nn/32,25,13/
!     do 1 i=1,32
!     do 1 j=1,25
!     do 1 k=1,13
!  1  data(i,j,k)=complex value
!     call fourt(data,nn,3,-1,1,work)
!
!     example 2.  one-dimensional forward transform of a real array of
!     length 64 in fortran ii,
!     dimension data(2,64)
!     do 2 i=1,64
!     data(1,i)=real part
!  2  data(2,i)=0.
!     call fourt(data,64,1,-1,0,0)
!
!     there are no error messages or error halts in this program.  the
!     program returns immediately if ndim or any nn(i) is less than one.
!
!     program by norman brenner from the basic program by charles
!     rader,  june 1967.  the idea for the digit reversal was
!     suggested by ralph alter.
!
!     this is the fastest and most versatile version of the fft known
!     to the author.  a program called four2 is available that also
!     performs the fast fourier transform and is written in usasi basic
!     fortran.  it is about one third as long and restricts the
!     dimensions of the input array (which must be complex) to be powers
!     of two.  another program, called four1, is one tenth as long and
!     runs two thirds as fast on a one-dimensional complex array whose
!     length is a power of two.
!
!     reference--
!     ieee audio transactions (june 1967), special issue on the fft.

      INTEGER, INTENT(IN) :: nn(3),ndim,isign,iform
      COMPLEX(DPC), INTENT(INOUT) :: data3(nn(1),nn(2),nn(3))
      REAL(EB), INTENT (INOUT) :: work(2*nn(1))
      
      INTEGER :: ifact(32),ntot,idim,np1,n,np2,m,ntwo,iif,idiv,iquot,irem,inon2,     &
                 icase,ifmin,i1rng,i,j,k,np2hf,i2,i1max,i1,i3,j3,nwork,ifp2,ifp1,    &
                 i2max,np1tw,ipar,k1,k2,mmax,lmax,l,kmin,kdif,kstep,k3,k4,np1hf,     &
                 j1min,j1,j2min,j2max,j2,j2rng,j1max,j3max,jmin,jmax,iconj,nhalf,    &
                 imin,imax,nprev=0,np0=0

      REAL(EB) :: data(2*nn(1)*nn(2)*nn(3)),tempr,tempi,u1r,u1i,u2r,u2i,u3r,         &
                  u3i,u4r,u4i,t2r,t2i,t3r,t3i,t4r,t4i,sumr,sumi,oldsr,oldsi,         &
                  difr,difi,theta,wr,wi,w2r,w2i,w3r,w3i,wstpr,wstpi,twowr
      
      REAL(EB), PARAMETER :: twopi=6.2831853071796, rthlf=0.70710678118655
      
      ! reshape data3 to 1D array
      data=0._EB
      n=1
      do k=1,nn(3)
        do j=1,nn(2)
          do i=1,nn(1)
            data(n)=real(data3(i,j,k),EB)
            data(n+1)=aimag(data3(i,j,k))
            n=n+2
          enddo
        enddo
      enddo
      
      if(ndim-1)920,1,1
1     ntot=2
      do 2 idim=1,ndim
      if(nn(idim))920,920,2
2     ntot=ntot*nn(idim)
!
!     main loop for each dimension
!
      np1=2
      do 910 idim=1,ndim
      n=nn(idim)
      np2=np1*n
      if(n-1)920,900,5
!
!     is n a power of two and if not, what are its factors
!
5     m=n
      ntwo=np1
      iif=1
      idiv=2
10    iquot=m/idiv
      irem=m-idiv*iquot
      if(iquot-idiv)50,11,11
11    if(irem)20,12,20
12    ntwo=ntwo+ntwo
      ifact(iif)=idiv
      iif=iif+1
      m=iquot
      go to 10
20    idiv=3
      inon2=iif
30    iquot=m/idiv
      irem=m-idiv*iquot
      if(iquot-idiv)60,31,31
31    if(irem)40,32,40
32    ifact(iif)=idiv
      iif=iif+1
      m=iquot
      go to 30
40    idiv=idiv+2
      go to 30
50    inon2=iif
      if(irem)60,51,60
51    ntwo=ntwo+ntwo
      go to 70
60    ifact(iif)=m
!
!     separate four cases--
!        1. complex transform or real transform for the 4th, 9th,etc.
!           dimensions.
!        2. real transform for the 2nd or 3rd dimension.  method--
!           transform half the data, supplying the other half by con-
!           jugate symmetry.
!        3. real transform for the 1st dimension, n odd.  method--
!           set the imaginary parts to zero.
!        4. real transform for the 1st dimension, n even.  method--
!           transform a complex array of length n/2 whose real parts
!           are the even numbered real values and whose imaginary parts
!           are the odd numbered real values.  separate and supply
!           the second half by conjugate symmetry.
!
70    icase=1
      ifmin=1
      i1rng=np1
      if(idim-4)71,100,100
71    if(iform)72,72,100
72    icase=2
      i1rng=np0*(1+nprev/2)
      if(idim-1)73,73,100
73    icase=3
      i1rng=np1
      if(ntwo-np1)100,100,74
74    icase=4
      ifmin=2
      ntwo=ntwo/2
      n=n/2
      np2=np2/2
      ntot=ntot/2
      i=1
      do 80 j=1,ntot
      data(j)=data(i)
80    i=i+2
!
!     shuffle data by bit reversal, since n=2**k.  as the shuffling
!     can be done by simple interchange, no working array is needed
!
100   if(ntwo-np2)200,110,110
110   np2hf=np2/2
      j=1
      do 150 i2=1,np2,np1
      if(j-i2)120,130,130
120   i1max=i2+np1-2
      do 125 i1=i2,i1max,2
      do 125 i3=i1,ntot,np2
      j3=j+i3-i2
      tempr=data(i3)
      tempi=data(i3+1)
      data(i3)=data(j3)
      data(i3+1)=data(j3+1)
      data(j3)=tempr
125   data(j3+1)=tempi
130   m=np2hf
140   if(j-m)150,150,145
145   j=j-m
      m=m/2
      if(m-np1)150,140,140
150   j=j+m
      go to 300
!
!     shuffle data by digit reversal for general n
!
200   nwork=2*n
      do 270 i1=1,np1,2
      do 270 i3=i1,ntot,np2
      j=i3
      do 260 i=1,nwork,2
      if(icase-3)210,220,210
210   work(i)=data(j)
      work(i+1)=data(j+1)
      go to 230
220   work(i)=data(j)
      work(i+1)=0._EB
230   ifp2=np2
      iif=ifmin
240   ifp1=ifp2/ifact(iif)
      j=j+ifp1
      if(j-i3-ifp2)260,250,250
250   j=j-ifp2
      ifp2=ifp1
      iif=iif+1
      if(ifp2-np1)260,260,240
260   continue
      i2max=i3+np2-np1
      i=1
      do 270 i2=i3,i2max,np1
      data(i2)=work(i)
      data(i2+1)=work(i+1)
270   i=i+2
!
!     main loop for factors of two.  perform fourier transforms of
!     length four, with one of length two if needed.  the twiddle factor
!     w=exp(isign*2*pi*sqrt(-1)*m/(4*mmax)).  check for w=isign*sqrt(-1)
!     and repeat for w=w*(1+isign*sqrt(-1))/sqrt(2).
!
300   if(ntwo-np1)600,600,305
305   np1tw=np1+np1
      ipar=ntwo/np1
310   if(ipar-2)350,330,320
320   ipar=ipar/4
      go to 310
330   do 340 i1=1,i1rng,2
      do 340 k1=i1,ntot,np1tw
      k2=k1+np1
      tempr=data(k2)
      tempi=data(k2+1)
      data(k2)=data(k1)-tempr
      data(k2+1)=data(k1+1)-tempi
      data(k1)=data(k1)+tempr
340   data(k1+1)=data(k1+1)+tempi
350   mmax=np1
360   if(mmax-ntwo/2)370,600,600
370   lmax=max0(np1tw,mmax/2)
      do 570 l=np1,lmax,np1tw
      m=l
      if(mmax-np1)420,420,380
380   theta=-twopi*REAL(l,EB)/REAL(4*mmax,EB)
      if(isign)400,390,390
390   theta=-theta
400   wr=cos(theta)
      wi=sin(theta)
410   w2r=wr*wr-wi*wi
      w2i=2._EB*wr*wi
      w3r=w2r*wr-w2i*wi
      w3i=w2r*wi+w2i*wr
420   do 530 i1=1,i1rng,2
      kmin=i1+ipar*m
      if(mmax-np1)430,430,440
430   kmin=i1
440   kdif=ipar*mmax
450   kstep=4*kdif
      if(kstep-ntwo)460,460,530
460   do 520 k1=kmin,ntot,kstep
      k2=k1+kdif
      k3=k2+kdif
      k4=k3+kdif
      if(mmax-np1)470,470,480
470   u1r=data(k1)+data(k2)
      u1i=data(k1+1)+data(k2+1)
      u2r=data(k3)+data(k4)
      u2i=data(k3+1)+data(k4+1)
      u3r=data(k1)-data(k2)
      u3i=data(k1+1)-data(k2+1)
      if(isign)471,472,472
471   u4r=data(k3+1)-data(k4+1)
      u4i=data(k4)-data(k3)
      go to 510
472   u4r=data(k4+1)-data(k3+1)
      u4i=data(k3)-data(k4)
      go to 510
480   t2r=w2r*data(k2)-w2i*data(k2+1)
      t2i=w2r*data(k2+1)+w2i*data(k2)
      t3r=wr*data(k3)-wi*data(k3+1)
      t3i=wr*data(k3+1)+wi*data(k3)
      t4r=w3r*data(k4)-w3i*data(k4+1)
      t4i=w3r*data(k4+1)+w3i*data(k4)
      u1r=data(k1)+t2r
      u1i=data(k1+1)+t2i
      u2r=t3r+t4r
      u2i=t3i+t4i
      u3r=data(k1)-t2r
      u3i=data(k1+1)-t2i
      if(isign)490,500,500
490   u4r=t3i-t4i
      u4i=t4r-t3r
      go to 510
500   u4r=t4i-t3i
      u4i=t3r-t4r
510   data(k1)=u1r+u2r
      data(k1+1)=u1i+u2i
      data(k2)=u3r+u4r
      data(k2+1)=u3i+u4i
      data(k3)=u1r-u2r
      data(k3+1)=u1i-u2i
      data(k4)=u3r-u4r
520   data(k4+1)=u3i-u4i
      kdif=kstep
      kmin=4*(kmin-i1)+i1
      go to 450
530   continue
      m=m+lmax
      if(m-mmax)540,540,570
540   if(isign)550,560,560
550   tempr=wr
      wr=(wr+wi)*rthlf
      wi=(wi-tempr)*rthlf
      go to 410
560   tempr=wr
      wr=(wr-wi)*rthlf
      wi=(tempr+wi)*rthlf
      go to 410
570   continue
      ipar=3-ipar
      mmax=mmax+mmax
      go to 360
!
!     main loop for factors not equal to two.  apply the twiddle factor
!     w=exp(isign*2*pi*sqrt(-1)*(j1-1)*(j2-j1)/(ifp1+ifp2)), then
!     perform a fourier transform of length ifact(iif), making use of
!     conjugate symmetries.
!
600   if(ntwo-np2)605,700,700
605   ifp1=ntwo
      iif=inon2
      np1hf=np1/2
610   ifp2=ifact(iif)*ifp1
      j1min=np1+1
      if(j1min-ifp1)615,615,640
615   do 635 j1=j1min,ifp1,np1
      theta=-twopi*REAL(j1-1,EB)/REAL(ifp2,EB)
      if(isign)625,620,620
620   theta=-theta
625   wstpr=cos(theta)
      wstpi=sin(theta)
      wr=wstpr
      wi=wstpi
      j2min=j1+ifp1
      j2max=j1+ifp2-ifp1
      do 635 j2=j2min,j2max,ifp1
      i1max=j2+i1rng-2
      do 630 i1=j2,i1max,2
      do 630 j3=i1,ntot,ifp2
      tempr=data(j3)
      data(j3)=data(j3)*wr-data(j3+1)*wi
630   data(j3+1)=tempr*wi+data(j3+1)*wr
      tempr=wr
      wr=wr*wstpr-wi*wstpi
635   wi=tempr*wstpi+wi*wstpr
640   theta=-twopi/REAL(ifact(iif),EB)
      if(isign)650,645,645
645   theta=-theta
650   wstpr=cos(theta)
      wstpi=sin(theta)
      j2rng=ifp1*(1+ifact(iif)/2)
      do 695 i1=1,i1rng,2
      do 695 i3=i1,ntot,np2
      j2max=i3+j2rng-ifp1
      do 690 j2=i3,j2max,ifp1
      j1max=j2+ifp1-np1
      do 680 j1=j2,j1max,np1
      j3max=j1+np2-ifp2
      do 680 j3=j1,j3max,ifp2
      jmin=j3-j2+i3
      jmax=jmin+ifp2-ifp1
      i=1+(j3-i3)/np1hf
      if(j2-i3)655,655,665
655   sumr=0._EB
      sumi=0._EB
      do 660 j=jmin,jmax,ifp1
659   sumr=sumr+data(j)
660   sumi=sumi+data(j+1)
      work(i)=sumr
      work(i+1)=sumi
      go to 680
665   iconj=1+(ifp2-2*j2+i3+j3)/np1hf
      j=jmax
      sumr=data(j)
      sumi=data(j+1)
      oldsr=0._EB
      oldsi=0._EB
      j=j-ifp1
670   tempr=sumr
      tempi=sumi
      sumr=twowr*sumr-oldsr+data(j)
      sumi=twowr*sumi-oldsi+data(j+1)
      oldsr=tempr
      oldsi=tempi
      j=j-ifp1
      if(j-jmin)675,675,670
675   tempr=wr*sumr-oldsr+data(j)
      tempi=wi*sumi
      work(i)=tempr-tempi
      work(iconj)=tempr+tempi
      tempr=wr*sumi-oldsi+data(j+1)
      tempi=wi*sumr
      work(i+1)=tempr+tempi
      work(iconj+1)=tempr-tempi
680   continue
      if(j2-i3)685,685,686
685   wr=wstpr
      wi=wstpi
      go to 690
686   tempr=wr
      wr=wr*wstpr-wi*wstpi
      wi=tempr*wstpi+wi*wstpr
690   twowr=wr+wr
      i=1
      i2max=i3+np2-np1
      do 695 i2=i3,i2max,np1
      data(i2)=work(i)
      data(i2+1)=work(i+1)
695   i=i+2
      iif=iif+1
      ifp1=ifp2
      if(ifp1-np2)610,700,700
!
!     complete a real transform in the 1st dimension, n even, by con-
!     jugate symmetries.
!
700   go to (900,800,900,701),icase
701   nhalf=n
      n=n+n
      theta=-twopi/REAL(n,EB)
      if(isign)703,702,702
702   theta=-theta
703   wstpr=cos(theta)
      wstpi=sin(theta)
      wr=wstpr
      wi=wstpi
      imin=3
      jmin=2*nhalf-1
      go to 725
710   j=jmin
      do 720 i=imin,ntot,np2
      sumr=(data(i)+data(j))/2._EB
      sumi=(data(i+1)+data(j+1))/2._EB
      difr=(data(i)-data(j))/2._EB
      difi=(data(i+1)-data(j+1))/2._EB
      tempr=wr*sumi+wi*difr
      tempi=wi*sumi-wr*difr
      data(i)=sumr+tempr
      data(i+1)=difi+tempi
      data(j)=sumr-tempr
      data(j+1)=-difi+tempi
720   j=j+np2
      imin=imin+2
      jmin=jmin-2
      tempr=wr
      wr=wr*wstpr-wi*wstpi
      wi=tempr*wstpi+wi*wstpr
725   if(imin-jmin)710,730,740
730   if(isign)731,740,740
731   do 735 i=imin,ntot,np2
735   data(i+1)=-data(i+1)
740   np2=np2+np2
      ntot=ntot+ntot
      j=ntot+1
      imax=ntot/2+1
745   imin=imax-2*nhalf
      i=imin
      go to 755
750   data(j)=data(i)
      data(j+1)=-data(i+1)
755   i=i+2
      j=j-2
      if(i-imax)750,760,760
760   data(j)=data(imin)-data(imin+1)
      data(j+1)=0._EB
      if(i-j)770,780,780
765   data(j)=data(i)
      data(j+1)=data(i+1)
770   i=i-2
      j=j-2
      if(i-imin)775,775,765
775   data(j)=data(imin)+data(imin+1)
      data(j+1)=0._EB
      imax=imin
      go to 745
780   data(1)=data(1)+data(2)
      data(2)=0._EB
      go to 900
!
!     complete a real transform for the 2nd or 3rd dimension by
!     conjugate symmetries.
!
800   if(i1rng-np1)805,900,900
805   do 860 i3=1,ntot,np2
      i2max=i3+np2-np1
      do 860 i2=i3,i2max,np1
      imin=i2+i1rng
      imax=i2+np1-2
      jmax=2*i3+np1-imin
      if(i2-i3)820,820,810
810   jmax=jmax+np2
820   if(idim-2)850,850,830 ! if (idim-2=1) goto 850; if (idim-2=2) goto 850; if (idim-2=3) goto 830
830   j=jmax+np0
      do 840 i=imin,imax,2
      data(i)=data(j)
      data(i+1)=-data(j+1)
840   j=j-2
850   j=jmax
      do 860 i=imin,imax,np0
      data(i)=data(j)
      data(i+1)=-data(j+1)
860   j=j-np0
!
!     end of loop on each dimension
!
900   np0=np1
      np1=np2
910   nprev=n

      ! reshape data back to 3D complex array
      
      !! for debug purposes (move to 920)
      !print *,size(data)
      !do i=1,size(data)
      !   print *,data(i)
      !enddo
      !stop
      
920   n=1
      do k=1,nn(3)
        do j=1,nn(2)
          do i=1,nn(1)
            data3(i,j,k)=cmplx(data(n),data(n+1),kind=DPC)
            n=n+2
          enddo
        enddo
      enddo                         
      return
      
END SUBROUTINE fourt_f90


SUBROUTINE ANALYTICAL_SOLUTION(NM)
IMPLICIT NONE
! Initialize flow variables with an analytical solution of the governing equations

INTEGER, INTENT(IN) :: NM
INTEGER :: I,J,K
REAL(EB) :: UU,WW

CALL POINT_TO_MESH(NM)

DO K=1,KBAR
   DO J=1,JBAR
      DO I=0,IBAR
         U(I,J,K) = 1._EB - 2._EB*COS(X(I))*SIN(ZC(K))
      ENDDO
   ENDDO
ENDDO
DO K=0,KBAR
   DO J=1,JBAR
      DO I=1,IBAR
         W(I,J,K) = 1._EB + 2._EB*SIN(XC(I))*COS(Z(K))
      ENDDO
   ENDDO
ENDDO
DO K=0,KBP1
   DO J=0,JBP1
      DO I=0,IBP1
         UU = 1._EB - 2._EB*COS(XC(I))*SIN(ZC(K))
         WW = 1._EB + 2._EB*SIN(XC(I))*COS(ZC(K))
         H(I,J,K) = -( COS(2._EB*XC(I)) + COS(2._EB*ZC(K)) ) + 0.5_EB*(UU**2+WW**2)
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE ANALYTICAL_SOLUTION


SUBROUTINE COMPRESSION_WAVE(NM,T,ITEST)
IMPLICIT NONE

INTEGER, INTENT(IN) :: NM,ITEST
REAL(EB), INTENT(IN) :: T
INTEGER :: I,J,K

CALL POINT_TO_MESH(NM)

SELECT CASE(ITEST)
   CASE(3) ! stationary compression wave
      DO K=1,KBAR
         DO J=1,JBAR
            DO I=0,IBAR
               U(I,J,K)  = 2._EB + SIN(X(I))
               US(I,J,K) = U(I,J,K)
            ENDDO
         ENDDO
      ENDDO
      DO K=1,KBAR
         DO J=0,JBAR
            DO I=1,IBAR
               V(I,J,K)  = 0._EB
               VS(I,J,K) = 0._EB
            ENDDO
         ENDDO
      ENDDO
      DO K=0,KBAR
         DO J=1,JBAR
            DO I=1,IBAR
               W(I,J,K)  = 3._EB + SIN(Z(K))
               WS(I,J,K) = W(I,J,K)
            ENDDO
         ENDDO
      ENDDO
   CASE(4) ! pulsating dilation
      DO K=1,KBAR
         DO J=1,JBAR
            DO I=0,IBAR
               U(I,J,K)  = 0.5_EB*SIN(X(I))*COS(T)
               US(I,J,K) = U(I,J,K)
            ENDDO
         ENDDO
      ENDDO
      DO K=1,KBAR
         DO J=0,JBAR
            DO I=1,IBAR
               V(I,J,K)  = 0._EB
               VS(I,J,K) = 0._EB
            ENDDO
         ENDDO
      ENDDO
      DO K=0,KBAR
         DO J=1,JBAR
            DO I=1,IBAR
               W(I,J,K)  = 0.5_EB*SIN(Z(K))*COS(T)
               WS(I,J,K) = W(I,J,K)
            ENDDO
         ENDDO
      ENDDO
END SELECT

DO K=1,KBAR
   DO J=1,JBAR
      DO I=1,IBAR
         D(I,J,K) = (U(I,J,K)-U(I-1,J,K))*RDX(I) + (W(I,J,K)-W(I,J,K-1))*RDZ(K)
         DS(I,J,K) = D(I,J,K)
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE COMPRESSION_WAVE


SUBROUTINE VARDEN_DYNSMAG(NM)
IMPLICIT NONE

!--------------------------------------------------------------
!     for all tensors, the indices are defined as follows...
!
!     |  11    12    13  |
!     |                  |
!     |  21    22    23  |
!     |                  |
!     |  31    32    33  |
!
!     I"name"I is the 'magnitude' of "name".
!--------------------------------------------------------------

INTEGER, INTENT(IN) :: NM

! Velocities relative to the p-cell center
REAL(EB) :: U_E,U_W,U_N,U_S,U_T,U_B
REAL(EB) :: V_E,V_W,V_N,V_S,V_T,V_B
REAL(EB) :: W_E,W_W,W_N,W_S,W_T,W_B
REAL(EB) :: DELTA,SKK
INTEGER :: I,J,K,N_LO(3),N_HI(3),ARRAY_LO(3),ARRAY_HI(3)

REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,UP,VP,WP,RHOP,RHOPHAT
REAL(EB), POINTER, DIMENSION(:,:,:) :: S11,S22,S33,S12,S13,S23,SS
REAL(EB), POINTER, DIMENSION(:,:,:) :: SHAT11,SHAT22,SHAT33,SHAT12,SHAT13,SHAT23,SSHAT
REAL(EB), POINTER, DIMENSION(:,:,:) :: BETA11,BETA22,BETA33,BETA12,BETA13,BETA23
REAL(EB), POINTER, DIMENSION(:,:,:) :: BETAHAT11,BETAHAT22,BETAHAT33,BETAHAT12,BETAHAT13,BETAHAT23
REAL(EB), POINTER, DIMENSION(:,:,:) :: M11,M22,M33,M12,M13,M23,MM,MMHAT
REAL(EB), POINTER, DIMENSION(:,:,:) :: L11,L22,L33,L12,L13,L23,ML,MLHAT

REAL(EB), PARAMETER :: ALPHA = 6.0_EB ! See Lund, 1997 CTR briefs.

! References:
!
! M. Germano, U. Piomelli, P. Moin, and W. Cabot.  A dynamic subgrid-scale eddy viscosity model.
! Phys. Fluids A, 3(7):1760-1765, 1991.
!
! M. Pino Martin, U. Piomelli, and G. Candler. Subgrid-scale models for compressible large-eddy
! simulation. Theoret. Comput. Fluid Dynamics, 13:361-376, 2000.
!
! P. Moin, K. Squires, W. Cabot, and S. Lee.  A dynamic subgrid-scale model for compressible
! turbulence and scalar transport. Phys. Fluids A, 3(11):2746-2757, 1991.
!
! T. S. Lund. On the use of discrete filters for large eddy simulation.  Center for Turbulence
! Research Annual Research Briefs, 1997.
!
! R. McDermott. Variable density formulation of the dynamic Smagorinsky model.
! http://randy.mcdermott.googlepages.com/dynsmag_comp.pdf

! *****************************************************************************
! CAUTION WHEN MODIFYING: The order in which the tensor components are computed
! is important because we overwrite pointers several times to conserve memory.
! *****************************************************************************

CALL POINT_TO_MESH(NM)

N_LO = 1
N_HI = (/IBAR,JBAR,KBAR/)

ARRAY_LO = 0
ARRAY_HI = (/IBP1,JBP1,KBP1/)

IF (PREDICTOR) THEN
   UU=>U
   VV=>V
   WW=>W
   RHOP=>RHO
ELSE
   UU=>US
   VV=>VS
   WW=>WS
   RHOP=>RHOS
ENDIF

UP => TURB_WORK1
VP => TURB_WORK2
WP => TURB_WORK3

S11 => WORK1
S22 => WORK2
S33 => WORK3
S12 => WORK4
S13 => WORK5
S23 => WORK6
SS  => WORK7


!$OMP PARALLEL DO COLLAPSE(3) PRIVATE(K,J,I,U_E,U_W,U_N,U_S,U_T,U_B,V_E,V_W,V_N,V_S,V_T,V_B,W_E,W_W,W_N,W_S,W_T,W_B,SKK)
DO K = N_LO(3),N_HI(3)
   DO J = N_LO(2),N_HI(2)
      DO I = N_LO(1),N_HI(1)

         U_E = UU(I,J,K)
         U_W = UU(I-1,J,K)
         U_N = 0.25_EB*( UU(I,J,K) + UU(I-1,J,K) + UU(I,J+1,K) + UU(I-1,J+1,K) )
         U_S = 0.25_EB*( UU(I,J,K) + UU(I-1,J,K) + UU(I,J-1,K) + UU(I-1,J-1,K) )
         U_T = 0.25_EB*( UU(I,J,K) + UU(I-1,J,K) + UU(I,J,K+1) + UU(I-1,J,K+1) )
         U_B = 0.25_EB*( UU(I,J,K) + UU(I-1,J,K) + UU(I,J,K-1) + UU(I-1,J,K-1) )

         V_N = VV(I,J,K)
         V_S = VV(I,J-1,K)
         V_E = 0.25_EB*( VV(I,J,K) + VV(I,J-1,K) + VV(I+1,J,K) + VV(I+1,J-1,K) )
         V_W = 0.25_EB*( VV(I,J,K) + VV(I,J-1,K) + VV(I-1,J,K) + VV(I-1,J-1,K) )
         V_T = 0.25_EB*( VV(I,J,K) + VV(I,J-1,K) + VV(I,J,K+1) + VV(I,J-1,K+1) )
         V_B = 0.25_EB*( VV(I,J,K) + VV(I,J-1,K) + VV(I,J,K-1) + VV(I,J-1,K-1) )

         W_T = WW(I,J,K)
         W_B = WW(I,J,K-1)
         W_E = 0.25_EB*( WW(I,J,K) + WW(I,J,K-1) + WW(I+1,J,K) + WW(I+1,J,K-1) ) 
         W_W = 0.25_EB*( WW(I,J,K) + WW(I,J,K-1) + WW(I-1,J,K) + WW(I-1,J,K-1) )
         W_N = 0.25_EB*( WW(I,J,K) + WW(I,J,K-1) + WW(I,J+1,K) + WW(I,J+1,K-1) )
         W_S = 0.25_EB*( WW(I,J,K) + WW(I,J,K-1) + WW(I,J-1,K) + WW(I,J-1,K-1) )

         UP(I,J,K) = 0.5_EB*(U_E + U_W)
         VP(I,J,K) = 0.5_EB*(V_N + V_S)
         WP(I,J,K) = 0.5_EB*(W_T + W_B)

         ! calculate the grid strain rate tensor
         
         S11(I,J,K) = (U_E - U_W)/DX(I)
         S22(I,J,K) = (V_N - V_S)/DY(J)
         S33(I,J,K) = (W_T - W_B)/DZ(K)
         SKK = S11(I,J,K) + S22(I,J,K) + S33(I,J,K)
         S11(I,J,K) = S11(I,J,K) - ONTH*SKK
         S22(I,J,K) = S22(I,J,K) - ONTH*SKK
         S33(I,J,K) = S33(I,J,K) - ONTH*SKK
         S12(I,J,K) = 0.5_EB*( (U_N - U_S)/DY(J) + (V_E - V_W)/DX(I) )
         S13(I,J,K) = 0.5_EB*( (U_T - U_B)/DZ(K) + (W_E - W_W)/DX(I) )
         S23(I,J,K) = 0.5_EB*( (V_T - V_B)/DZ(K) + (W_N - W_S)/DY(J) )
         
         ! calculate magnitude of the grid strain rate

         SS(I,J,K) = SQRT(2._EB*(S11(I,J,K)*S11(I,J,K) + &
                                 S22(I,J,K)*S22(I,J,K) + &
                                 S33(I,J,K)*S33(I,J,K) + &
                          2._EB*(S12(I,J,K)*S12(I,J,K) + &
                                 S13(I,J,K)*S13(I,J,K) + &
                                 S23(I,J,K)*S23(I,J,K)) ) )

      ENDDO
   ENDDO
ENDDO
!$OMP END PARALLEL DO

! test filter the strain rate

SHAT11 => TURB_WORK4
SHAT22 => TURB_WORK5
SHAT33 => TURB_WORK6
SHAT12 => TURB_WORK7
SHAT13 => TURB_WORK8
SHAT23 => TURB_WORK9
SSHAT  => TURB_WORK10

CALL TEST_FILTER(SHAT11,S11,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(SHAT22,S22,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(SHAT33,S33,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(SHAT12,S12,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(SHAT13,S13,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(SHAT23,S23,N_LO,N_HI,ARRAY_LO,ARRAY_HI)

! calculate magnitude of test filtered strain rate

!$OMP PARALLEL DO COLLAPSE(3) PRIVATE(K,J,I)
DO K = N_LO(3),N_HI(3)
   DO J = N_LO(2),N_HI(2)
      DO I = N_LO(1),N_HI(1)
      
         SSHAT(I,J,K) = SQRT(2._EB*(SHAT11(I,J,K)*SHAT11(I,J,K) + &
                                    SHAT22(I,J,K)*SHAT22(I,J,K) + &
                                    SHAT33(I,J,K)*SHAT33(I,J,K) + &
                             2._EB*(SHAT12(I,J,K)*SHAT12(I,J,K) + &
                                    SHAT13(I,J,K)*SHAT13(I,J,K) + &
                                    SHAT23(I,J,K)*SHAT23(I,J,K)) ) )

      ENDDO
   ENDDO
ENDDO
!$OMP END PARALLEL DO

! calculate the grid filtered stress tensor, beta

BETA11 => WORK1
BETA22 => WORK2
BETA33 => WORK3
BETA12 => WORK4
BETA13 => WORK5
BETA23 => WORK6

!$OMP PARALLEL WORKSHARE
BETA11 = RHOP*SS*S11
BETA22 = RHOP*SS*S22
BETA33 = RHOP*SS*S33
BETA12 = RHOP*SS*S12
BETA13 = RHOP*SS*S13
BETA23 = RHOP*SS*S23
!$OMP END PARALLEL WORKSHARE

! test filter the grid filtered stress tensor

BETAHAT11 => WORK1
BETAHAT22 => WORK2
BETAHAT33 => WORK3
BETAHAT12 => WORK4
BETAHAT13 => WORK5
BETAHAT23 => WORK6

CALL TEST_FILTER(BETAHAT11,BETA11,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(BETAHAT22,BETA22,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(BETAHAT33,BETA33,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(BETAHAT12,BETA12,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(BETAHAT13,BETA13,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(BETAHAT23,BETA23,N_LO,N_HI,ARRAY_LO,ARRAY_HI)

! test filter the density

RHOPHAT => WORK7
CALL TEST_FILTER(RHOPHAT,RHOP,N_LO,N_HI,ARRAY_LO,ARRAY_HI)

! calculate the Mij tensor

M11 => WORK1
M22 => WORK2
M33 => WORK3
M12 => WORK4
M13 => WORK5
M23 => WORK6

!$OMP PARALLEL DO COLLAPSE(3) PRIVATE(K,J,I)
DO K = N_LO(3),N_HI(3)
   DO J = N_LO(2),N_HI(2)
      DO I = N_LO(1),N_HI(1)
      
         M11(I,J,K) = 2._EB*(BETAHAT11(I,J,K) - ALPHA*RHOPHAT(I,J,K)*SSHAT(I,J,K)*SHAT11(I,J,K))
         M22(I,J,K) = 2._EB*(BETAHAT22(I,J,K) - ALPHA*RHOPHAT(I,J,K)*SSHAT(I,J,K)*SHAT22(I,J,K))
         M33(I,J,K) = 2._EB*(BETAHAT33(I,J,K) - ALPHA*RHOPHAT(I,J,K)*SSHAT(I,J,K)*SHAT33(I,J,K))
         M12(I,J,K) = 2._EB*(BETAHAT12(I,J,K) - ALPHA*RHOPHAT(I,J,K)*SSHAT(I,J,K)*SHAT12(I,J,K))
         M13(I,J,K) = 2._EB*(BETAHAT13(I,J,K) - ALPHA*RHOPHAT(I,J,K)*SSHAT(I,J,K)*SHAT13(I,J,K))
         M23(I,J,K) = 2._EB*(BETAHAT23(I,J,K) - ALPHA*RHOPHAT(I,J,K)*SSHAT(I,J,K)*SHAT23(I,J,K))
         
      ENDDO
   ENDDO
ENDDO
!$OMP END PARALLEL DO

! calculate the Leonard term, Lij

L11 => TURB_WORK4
L22 => TURB_WORK5
L33 => TURB_WORK6
L12 => TURB_WORK7
L13 => TURB_WORK8
L23 => TURB_WORK9

CALL CALC_VARDEN_LEONARD_TERM(NM)

! calculate Mij*Lij & Mij*Mij

MM    => TURB_WORK1
MMHAT => TURB_WORK1

ML    => TURB_WORK2
MLHAT => TURB_WORK2

!$OMP PARALLEL DO COLLAPSE(3) PRIVATE(K,J,I)
DO K = N_LO(3),N_HI(3)
   DO J = N_LO(2),N_HI(2)
      DO I = N_LO(1),N_HI(1)
      
         ML(I,J,K) = M11(I,J,K)*L11(I,J,K) + M22(I,J,K)*L22(I,J,K) + M33(I,J,K)*L33(I,J,K) + &
              2._EB*(M12(I,J,K)*L12(I,J,K) + M13(I,J,K)*L13(I,J,K) + M23(I,J,K)*L23(I,J,K))
       
         MM(I,J,K) = M11(I,J,K)*M11(I,J,K) + M22(I,J,K)*M22(I,J,K) + M33(I,J,K)*M33(I,J,K) + &
              2._EB*(M12(I,J,K)*M12(I,J,K) + M13(I,J,K)*M13(I,J,K) + M23(I,J,K)*M23(I,J,K))
              
      ENDDO
   ENDDO
ENDDO
!$OMP END PARALLEL DO

! do some smoothing

CALL TEST_FILTER(MLHAT,ML,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(MMHAT,MM,N_LO,N_HI,ARRAY_LO,ARRAY_HI)

!$OMP PARALLEL DO COLLAPSE(3) PRIVATE(K,J,I,DELTA)
DO K = N_LO(3),N_HI(3)
   DO J = N_LO(2),N_HI(2)
      DO I = N_LO(1),N_HI(1)

         ! calculate the local Smagorinsky coefficient

         ! perform "clipping" in case MLij is negative...
         IF (MLHAT(I,J,K) < 0._EB) MLHAT(I,J,K) = 0._EB

         ! calculate the effective viscosity

         ! handle the case where we divide by zero, note MMHAT is positive semi-definite
         IF (MMHAT(I,J,K) == 0._EB) THEN
            C_DYNSMAG(I,J,K) = 0._EB
         ELSE
            ! filter width
            IF (TWO_D) THEN
               DELTA = (DX(I)*DZ(K))**0.5_EB
            ELSE
               DELTA = (DX(I)*DY(J)*DZ(K))**ONTH
            ENDIF
            C_DYNSMAG(I,J,K) = SQRT(MLHAT(I,J,K)/MMHAT(I,J,K))/DELTA
         ENDIF
         
         ! clip max value of CS, note that CS*DELTA is the "mixing length", so DELTA
         ! is a reasonable upper bound
         C_DYNSMAG(I,J,K) = MIN(C_DYNSMAG(I,J,K),1._EB)
				
      END DO
   END DO
END DO
!$OMP END PARALLEL DO

END SUBROUTINE VARDEN_DYNSMAG


SUBROUTINE CALC_VARDEN_LEONARD_TERM(NM)
IMPLICIT NONE

INTEGER, INTENT(IN) :: NM

REAL(EB), POINTER, DIMENSION(:,:,:) :: L11,L22,L33,L12,L13,L23
REAL(EB), POINTER, DIMENSION(:,:,:) :: RHOP,RHOPHAT
REAL(EB), POINTER, DIMENSION(:,:,:) :: UP,VP,WP
REAL(EB), POINTER, DIMENSION(:,:,:) :: RUU,RVV,RWW,RUV,RUW,RVW
REAL(EB), POINTER, DIMENSION(:,:,:) :: RU,RV,RW
REAL(EB), POINTER, DIMENSION(:,:,:) :: RUU_HAT,RVV_HAT,RWW_HAT,RUV_HAT,RUW_HAT,RVW_HAT
REAL(EB), POINTER, DIMENSION(:,:,:) :: RU_HAT,RV_HAT,RW_HAT

TYPE(MESH_TYPE), POINTER :: M
INTEGER :: I,J,K,N_LO(3),N_HI(3),ARRAY_LO(3),ARRAY_HI(3)

! *****************************************************************************
! CAUTION WHEN MODIFYING: The order in which the tensor components are computed
! is important because we overwrite pointers several times to conserve memory.
! *****************************************************************************

M => MESHES(NM)

IF (PREDICTOR) THEN
   RHOP=>RHO
ELSE
   RHOP=>RHOS
ENDIF
RHOPHAT => M%WORK7

N_LO = 1
N_HI = (/M%IBAR,M%JBAR,M%KBAR/)

ARRAY_LO = 0
ARRAY_HI = (/M%IBP1,M%JBP1,M%KBP1/)

! Compute rho*UiUj

UP => M%TURB_WORK1 ! will be overwritten by RU
VP => M%TURB_WORK2
WP => M%TURB_WORK3

RUU => M%TURB_WORK4 ! will be overwritten by RUU_HAT
RVV => M%TURB_WORK5
RWW => M%TURB_WORK6
RUV => M%TURB_WORK7
RUW => M%TURB_WORK8
RVW => M%TURB_WORK9

!$OMP PARALLEL DO COLLAPSE(3) PRIVATE(K,J,I)
DO K=N_LO(3),N_HI(3)
   DO J=N_LO(2),N_HI(2)
      DO I=N_LO(1),N_HI(1)

         RUU(I,J,K) = RHOP(I,J,K)*UP(I,J,K)*UP(I,J,K)
         RVV(I,J,K) = RHOP(I,J,K)*VP(I,J,K)*VP(I,J,K)
         RWW(I,J,K) = RHOP(I,J,K)*WP(I,J,K)*WP(I,J,K)
         RUV(I,J,K) = RHOP(I,J,K)*UP(I,J,K)*VP(I,J,K)
         RUW(I,J,K) = RHOP(I,J,K)*UP(I,J,K)*WP(I,J,K)
         RVW(I,J,K) = RHOP(I,J,K)*VP(I,J,K)*WP(I,J,K)

      ENDDO
   ENDDO
ENDDO
!$OMP END PARALLEL DO

! Test filter rho*UiUj

RUU_HAT => M%TURB_WORK4 ! will be overwritten by Lij
RVV_HAT => M%TURB_WORK5
RWW_HAT => M%TURB_WORK6
RUV_HAT => M%TURB_WORK7
RUW_HAT => M%TURB_WORK8
RVW_HAT => M%TURB_WORK9

CALL TEST_FILTER(RUU_HAT,RUU,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(RVV_HAT,RVV,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(RWW_HAT,RWW,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(RUV_HAT,RUV,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(RUW_HAT,RUW,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(RVW_HAT,RVW,N_LO,N_HI,ARRAY_LO,ARRAY_HI)

! Compute rho*Ui

RU => M%TURB_WORK1 ! will be overwritten by RU_HAT
RV => M%TURB_WORK2
RW => M%TURB_WORK3

!$OMP PARALLEL WORKSHARE
RU = RHOP*UP
RV = RHOP*VP
RW = RHOP*WP
!$OMP END PARALLEL WORKSHARE

! Test filter rho*Ui

RU_HAT => M%TURB_WORK1
RV_HAT => M%TURB_WORK2
RW_HAT => M%TURB_WORK3

CALL TEST_FILTER(RU_HAT,RU,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(RV_HAT,RV,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(RW_HAT,RW,N_LO,N_HI,ARRAY_LO,ARRAY_HI)

! Compute variable density Leonard stress

L11 => M%TURB_WORK4
L22 => M%TURB_WORK5
L33 => M%TURB_WORK6
L12 => M%TURB_WORK7
L13 => M%TURB_WORK8
L23 => M%TURB_WORK9

!$OMP PARALLEL DO COLLAPSE(3) PRIVATE(K,J,I)
DO K=N_LO(3),N_HI(3)
   DO J=N_LO(2),N_HI(2)
      DO I=N_LO(1),N_HI(1)
         IF (RHOPHAT(I,J,K)>0._EB) THEN
            L11(I,J,K) = RUU_HAT(I,J,K) - RU_HAT(I,J,K)*RU_HAT(I,J,K)/RHOPHAT(I,J,K)
            L22(I,J,K) = RVV_HAT(I,J,K) - RV_HAT(I,J,K)*RV_HAT(I,J,K)/RHOPHAT(I,J,K)
            L33(I,J,K) = RWW_HAT(I,J,K) - RW_HAT(I,J,K)*RW_HAT(I,J,K)/RHOPHAT(I,J,K)
            L12(I,J,K) = RUV_HAT(I,J,K) - RU_HAT(I,J,K)*RV_HAT(I,J,K)/RHOPHAT(I,J,K)
            L13(I,J,K) = RUW_HAT(I,J,K) - RU_HAT(I,J,K)*RW_HAT(I,J,K)/RHOPHAT(I,J,K)
            L23(I,J,K) = RVW_HAT(I,J,K) - RV_HAT(I,J,K)*RW_HAT(I,J,K)/RHOPHAT(I,J,K)
         ELSE
            L11(I,J,K) = 0._EB
            L22(I,J,K) = 0._EB
            L33(I,J,K) = 0._EB
            L12(I,J,K) = 0._EB
            L13(I,J,K) = 0._EB
            L23(I,J,K) = 0._EB
         ENDIF
      ENDDO
   ENDDO
ENDDO
!$OMP END PARALLEL DO

END SUBROUTINE CALC_VARDEN_LEONARD_TERM


SUBROUTINE TEST_FILTER(PHIBAR,PHI,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
IMPLICIT NONE

INTEGER, INTENT(IN) :: N_LO(3),N_HI(3),ARRAY_LO(3),ARRAY_HI(3)
REAL(EB), INTENT(IN) :: PHI(ARRAY_LO(1):ARRAY_HI(1),ARRAY_LO(2):ARRAY_HI(2),ARRAY_LO(3):ARRAY_HI(3))
REAL(EB), INTENT(OUT) :: PHIBAR(ARRAY_LO(1):ARRAY_HI(1),ARRAY_LO(2):ARRAY_HI(2),ARRAY_LO(3):ARRAY_HI(3))
REAL(EB), POINTER, DIMENSION(:) :: PHI1,PHI2
INTEGER I,J,K

PHI1 => TURB_WORK11
PHI2 => TURB_WORK12

PHIBAR = PHI

! filter in x:
DO K = N_LO(3),N_HI(3)
   DO J = N_LO(2),N_HI(2)
      PHI1(N_LO(1):N_HI(1)) = PHIBAR(N_LO(1):N_HI(1),J,K)
      CALL TOPHAT_FILTER_1D(PHI2(N_LO(1):N_HI(1)),PHI1(N_LO(1):N_HI(1)),N_LO(1),N_HI(1))
      PHIBAR(N_LO(1):N_HI(1),J,K) = PHI2(N_LO(1):N_HI(1))
   ENDDO
ENDDO

IF (.NOT.TWO_D) THEN
   ! filter in y:
   DO K = N_LO(3),N_HI(3)
      DO I = N_LO(1),N_HI(1)
         PHI1(N_LO(2):N_HI(2)) = PHIBAR(I,N_LO(2):N_HI(2),K)
         CALL TOPHAT_FILTER_1D(PHI2(N_LO(2):N_HI(2)),PHI1(N_LO(2):N_HI(2)),N_LO(2),N_HI(2))
         PHIBAR(I,N_LO(2):N_HI(2),K) = PHI2(N_LO(2):N_HI(2))
      ENDDO
   ENDDO
ENDIF

! filter in z:
DO J = N_LO(2),N_HI(2)
   DO I = N_LO(1),N_HI(1)
      PHI1(N_LO(3):N_HI(3)) = PHIBAR(I,J,N_LO(3):N_HI(3))
      CALL TOPHAT_FILTER_1D(PHI2(N_LO(3):N_HI(3)),PHI1(N_LO(3):N_HI(3)),N_LO(3),N_HI(3))
      PHIBAR(I,J,N_LO(3):N_HI(3)) = PHI2(N_LO(3):N_HI(3))
   ENDDO
ENDDO

END SUBROUTINE TEST_FILTER


SUBROUTINE TOPHAT_FILTER_1D(UBAR,U,N_LO,N_HI)
IMPLICIT NONE

INTEGER, INTENT(IN) :: N_LO,N_HI
REAL(EB), INTENT(IN) :: U(N_LO:N_HI)
REAL(EB), INTENT(OUT) :: UBAR(N_LO:N_HI)
INTEGER :: J
REAL(EB), POINTER, DIMENSION(:) :: UU
REAL(EB):: W(-1:1)

UU => WORK
UU(N_LO:N_HI) = U

W(-1:1) = (/0.25_EB,0.5_EB,0.25_EB/)   ! trapezoid rule
!W(-1:1) = (/ONSI,TWTH,ONSI/)           ! Simpson's rule

! Filter the u field to obtain ubar
DO J=N_LO+1,N_HI-1
   UBAR(J) = W(-1)*UU(J-1) + W(0)*UU(J) + W(1)*UU(J+1)
ENDDO
! set boundary values (not ideal, but fast and simple)
UBAR(N_LO) = UBAR(N_LO+1)
UBAR(N_HI) = UBAR(N_HI-1)

END SUBROUTINE TOPHAT_FILTER_1D


SUBROUTINE TURBULENT_KINETIC_ENERGY(NM)
IMPLICIT NONE

INTEGER, INTENT(IN) :: NM
REAL(EB) :: K_TOT,K_LES,K_SGS
INTEGER :: I,J,K,N_LO(3),N_HI(3),ARRAY_LO(3),ARRAY_HI(3)
REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,UP,VP,WP,UP_HAT,VP_HAT,WP_HAT

CALL POINT_TO_MESH(NM)

N_LO = 1
N_HI = (/IBAR,JBAR,KBAR/)

ARRAY_LO = 0
ARRAY_HI = (/IBP1,JBP1,KBP1/)

IF (PREDICTOR) THEN
   UU=>US
   VV=>VS
   WW=>WS
ELSE
   UU=>U
   VV=>V
   WW=>W
ENDIF

! Velocities relative to the p-cell center
UP => WORK1
VP => WORK2
WP => WORK3

DO K = N_LO(3),N_HI(3)
   DO J = N_LO(2),N_HI(2)
      DO I = N_LO(1),N_HI(1)
         UP(I,J,K) = 0.5_EB*(UU(I,J,K) + UU(I-1,J,K))
         VP(I,J,K) = 0.5_EB*(VV(I,J,K) + VV(I,J-1,K))
         WP(I,J,K) = 0.5_EB*(WW(I,J,K) + WW(I,J,K-1))
      ENDDO
   ENDDO
ENDDO

UP_HAT => WORK4
VP_HAT => WORK5
WP_HAT => WORK6

CALL TEST_FILTER(UP_HAT,UP,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(VP_HAT,VP,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
CALL TEST_FILTER(WP_HAT,WP,N_LO,N_HI,ARRAY_LO,ARRAY_HI)

DO K = N_LO(3),N_HI(3)
   DO J = N_LO(2),N_HI(2)
      DO I = N_LO(1),N_HI(1)
      
         K_LES = UP(I,J,K)**2 + VP(I,J,K)**2 + WP(I,J,K)**2
         K_SGS = (UP(I,J,K)-UP_HAT(I,J,K))**2 + (VP(I,J,K)-VP_HAT(I,J,K))**2 + (WP(I,J,K)-WP_HAT(I,J,K))**2
         K_TOT = K_LES + K_SGS
         
         IF (K_TOT>0._EB) THEN
            MTR(I,J,K) = K_SGS/K_TOT
         ELSE
            MTR(I,J,K) = 0._EB
         ENDIF
         
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE TURBULENT_KINETIC_ENERGY


SUBROUTINE SUBGRID_SCALAR_VARIANCE(NM)
IMPLICIT NONE

INTEGER, INTENT(IN) :: NM
INTEGER :: I,J,K,N_LO(3),N_HI(3),ARRAY_LO(3),ARRAY_HI(3)
REAL(EB) :: T_TOT,T_LES,T_SGS
REAL(EB), POINTER, DIMENSION(:,:,:,:) :: YYP
REAL(EB), POINTER, DIMENSION(:,:,:) :: ZP,ZP_HAT

CALL POINT_TO_MESH(NM)

N_LO = 1
N_HI = (/IBAR,JBAR,KBAR/)

ARRAY_LO = 0
ARRAY_HI = (/IBP1,JBP1,KBP1/)

IF (PREDICTOR) THEN
   YYP=>YYS
ELSE
   YYP=>YY
ENDIF
ZP => WORK1
ZP_HAT => WORK2

DO K = ARRAY_LO(3),ARRAY_HI(3)
   DO J = ARRAY_LO(2),ARRAY_HI(2)
      DO I = ARRAY_LO(1),ARRAY_HI(1)
         ZP(I,J,K) = YYP(I,J,K,1) ! this is a hack for now because mixture fraction is not general
      ENDDO
   ENDDO
ENDDO
CALL TEST_FILTER(ZP_HAT,ZP,N_LO,N_HI,ARRAY_LO,ARRAY_HI)
   
DO K = N_LO(3),N_HI(3)
   DO J = N_LO(2),N_HI(2)
      DO I = N_LO(1),N_HI(1)
      
         T_LES = ZP(I,J,K)**2 - ZP_HAT(I,J,K)**2
         T_SGS = (ZP(I,J,K)-ZP_HAT(I,J,K))**2
         T_TOT = T_LES + T_SGS
         
         IF (T_TOT>0._EB) THEN
            MSR(I,J,K) = T_SGS/T_TOT
         ELSE
            MSR(I,J,K) = 0._EB
         ENDIF
         
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE SUBGRID_SCALAR_VARIANCE


SUBROUTINE WAVELET_ERROR_ESTIMATOR(NM)
IMPLICIT NONE

INTEGER, INTENT(IN) :: NM
INTEGER :: I,J,K
REAL(EB) :: SS(4)
REAL(EB), POINTER, DIMENSION(:,:,:,:) :: YYP
REAL(EB), POINTER, DIMENSION(:,:,:) :: ZP

CALL POINT_TO_MESH(NM)

IF (PREDICTOR) THEN
   YYP=>YYS
ELSE
   YYP=>YY
ENDIF
ZP => WORK1

DO K = 0,KBP1
   DO J = 0,JBP1
      DO I = 0,IBP1
         ZP(I,J,K) = YYP(I,J,K,1)
      ENDDO
   ENDDO
ENDDO

WEM = 0._EB
DO K = 1,KBAR
   DO J = 1,JBAR
      DO I = 2,IBAR
          SS = ZP(I-2:I+1,J,K)
          WEM(I,J,K) = WAVELET_ERROR(SS)
      ENDDO
   ENDDO
ENDDO
IF (.NOT.TWO_D) THEN
   DO K = 1,KBAR
      DO J = 2,JBAR
         DO I = 1,IBAR
            SS = ZP(I,J-2:J+1,K)
            WEM(I,J,K) = MAX(WEM(I,J,K),WAVELET_ERROR(SS))
         ENDDO
      ENDDO
   ENDDO
ENDIF
DO K = 2,KBAR
   DO J = 1,JBAR
      DO I = 1,IBAR
         SS = ZP(I,J,K-2:K+1)
         WEM(I,J,K) = MAX(WEM(I,J,K),WAVELET_ERROR(SS))
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE WAVELET_ERROR_ESTIMATOR


REAL(EB) FUNCTION WAVELET_ERROR(S)
IMPLICIT NONE

INTEGER, PARAMETER :: M=2 ! only need two level transform, but could be generalized
REAL(EB), INTENT(IN) :: S(2*M)
REAL(EB) :: A(M,M)=0._EB,C(M,M)=0._EB,C1,C2
INTEGER :: I,J,K,N

! Comments: This function generates a normalized error measure WAVELET_ERROR based on coefficients
! from a simple Haar wavelet transform.  The function requires the input of 5 scalar values.  The
! error is estimated at the point of the median value S(3) based on the 4 segment averages...
!
!    S(1)    S(2)    S(3)    S(4)
! |   o   |   o   |   o   |   o   |
!                     E

! discrete Haar wavelet transform
N=M
DO I=1,M
   DO J=1,N
      K=2*J-1
      IF (I==1) THEN
         A(I,J) = 0.5_EB*(S(K)+S(K+1))
         C(I,J) = 0.5_EB*(S(K)-S(K+1))
      ELSE
         A(I,J) = 0.5_EB*(A(I-1,K)+A(I-1,K+1))
         C(I,J) = 0.5_EB*(A(I-1,K)-A(I-1,K+1))
      ENDIF
   ENDDO
   N=N/2;
ENDDO

C1 = SUM(C(1,:))
C2 = SUM(C(2,:))

WAVELET_ERROR = ABS(C1-C2)/(ABS(A(2,1))+1.E-6_EB)

END FUNCTION WAVELET_ERROR


SUBROUTINE WERNER_WENGLE_WALL_MODEL(SF,U_TAU,U1,NU,DZ,ROUGHNESS)
IMPLICIT NONE

REAL(EB), INTENT(OUT) :: SF
REAL(EB), INTENT(IN) :: U1,NU,DZ,ROUGHNESS

REAL(EB), PARAMETER :: A=8.3_EB,B=1._EB/7._EB
REAL(EB), PARAMETER :: Z_PLUS_TURBULENT = 11.81_EB
REAL(EB), PARAMETER :: ALPHA=7.202125273562269_EB !! ALPHA=(1._EB-B)/2._EB*A**((1._EB+B)/(1._EB-B))
REAL(EB), PARAMETER :: BETA=1._EB+B
REAL(EB), PARAMETER :: ETA=(1._EB+B)/A
REAL(EB), PARAMETER :: GAMMA=2._EB/(1._EB+B)
REAL(EB), PARAMETER :: RKAPPA=2.44_EB ! 1./von Karman constant
REAL(EB), PARAMETER :: BTILDE=7.44_EB ! see Pope p. 297 (constant has been modified)

REAL(EB) :: U_TAU,TAU_W,NU_OVER_DZ,Z_PLUS,TAU_ROUGH

! References (for smooth walls):
!
! Werner, H., Wengle, H. (1991) Large-eddy simulation of turbulent flow over
! and around a cube in a plate channel. 8th Symposium on Turbulent Shear
! Flows, Munich, Germany.
!
! Pierre Sagaut. Large Eddy Simulation for Incompressible Flows: An Introduction.
! Springer, 2001.
!
! Temmerman, L., Leschziner, M.A., Mellen, C.P., and Frohlich, J. (2003)
! Investigation of wall-function approximations and subgrid-scale models in
! Large Eddy Simulation of separated flow in a channel with streamwise
! periodic constrictions. International Journal of Heat and Fluid Flow,
! Vol. 24, No. 2, pp. 157-180.
!
! Breuer, M., Kniazev, B., and Abel, M. (2007) Development of wall models
! for LES of separated flows using statistical evaluations. Computers and
! Fluids, Vol. 36, pp. 817-837.
!
! McDermott, R. (2009) FDS Wall Flows, Part I: Straight Channels, NIST Technical Note.
!
! References (for rough surfaces):
!
! S. B. Pope (2000) Turbulent Flows, Cambridge.
!
! Moeng, C.-H. (1984) A Large-Eddy Simulation Model for the Study of Planetary
! Boundary-Layer Turbulence. Journal of the Atmospheric Sciences, Vol. 41, No. 13,
! pp. 2052-2062.
!
! Stoll, R., Porte-Agel, F. (2008) Large-Eddy Simulation of the Stable Atmospheric
! Boundary Layer using Dynamic Models with Different Averaging Schemes. Boundary-Layer
! Meteorology, 126:1-28.
!
! Comments:
!
! The slip factor (SF) is based on the following approximation to the wall stress
! (note that u0 is the ghost cell value of the streamwise velocity component and
! z is the wall-normal direction):
! tau_w = mu*(u1-u0)/dz = mu*(u1-SF*u1)/dz = mu*u1/dz*(1-SF)
! note that tau_w/rho = nu*u1/dz*(1-SF)

TAU_ROUGH = 0._EB
IF (ROUGHNESS>0._EB) THEN
   ! Pope (2000)
   TAU_ROUGH = ( U1/(RKAPPA*LOG(0.5_EB*DZ/ROUGHNESS)+BTILDE) )**2 ! actually tau_w/rho
ENDIF
! Werner-Wengle
NU_OVER_DZ = NU/DZ
TAU_W = (ALPHA*(NU_OVER_DZ)**BETA + ETA*(NU_OVER_DZ)**B*ABS(U1))**GAMMA ! actually tau_w/rho
TAU_W = MAX(TAU_W,TAU_ROUGH)
U_TAU = SQRT(TAU_W)
Z_PLUS = DZ/(NU/(U_TAU+1.E-10_EB))
IF (Z_PLUS>Z_PLUS_TURBULENT) THEN
   SF = 1._EB-TAU_W/(NU/DZ*ABS(U1)) ! log layer
ELSE
   SF = -1._EB ! viscous sublayer
ENDIF

!! check values...
!IF (Z_PLUS>Z_PLUS_TURBULENT) THEN
!   print *,'A = ',A
!   print *,'B = ',B
!   print *,'ALPHA = ',ALPHA
!   print *,'BETA = ',BETA
!   print *,'ETA = ',ETA
!   print *,'GAMMA = ',GAMMA
!   print *,'U1 = ',U1
!   print *,'NU/DZ = ',NU_OVER_DZ
!   print *,'TAU_W/RHO = ',TAU_W
!   print *,'Z_PLUS = ',Z_PLUS
!   print *,'SF = ',SF
!   print *
!ENDIF

END SUBROUTINE WERNER_WENGLE_WALL_MODEL


SUBROUTINE SURFACE_HEAT_FLUX_MODEL(H,U_TAU,DZ,ROUGHNESS,IOR,RHO,CP)

REAL(EB), INTENT(OUT) :: H ! heat transfer coefficient
REAL(EB), INTENT(IN) :: U_TAU,DZ,ROUGHNESS,RHO,CP
INTEGER, INTENT(IN) :: IOR
REAL(EB), PARAMETER :: KAPPA=0.41_EB ! von Karman constant
REAL(EB) :: PSI,MOL,Z0

! References:
!
! Stoll, R., Porte-Agel, F. (2008) Large-Eddy Simulation of the Stable Atmospheric
! Boundary Layer using Dynamic Models with Different Averaging Schemes. Boundary-Layer
! Meteorology, 126:1-28.

PSI = 0._EB
MOL = 0._EB
Z0 = MAX(ROUGHNESS,1.E-6_EB)

! atmospheric stability correction (use later)
IF (IOR==3) THEN
   MOL = 0._EB !! -U_TAU**3*THETA/(KAPPA*GRAV*HEAT_FLUX)
   PSI = 0._EB !! -7.8_EB*0.5*DZ/MOL
ENDIF

H = RHO*CP*U_TAU*KAPPA/(LOG(0.5_EB*DZ/Z0)-PSI)

END SUBROUTINE SURFACE_HEAT_FLUX_MODEL


SUBROUTINE PERM3D(UP,VP,WP,Q,R,S,NC,UBAR,D1,D2,DIV)
	IMPLICIT NONE
	
	INTEGER, INTENT(IN) :: NC					! NUMBER OF PARTICLES IN THE CELL
	REAL, INTENT(IN) :: Q(NC),R(NC),S(NC)		! (Q,R,S) NORMALIZED POSITION OF PARTICLES
	REAL, INTENT(IN) :: UBAR(12),D1(12),D2(12)	! PARABOLIC EDGE PARAMETERS
	REAL, INTENT(OUT) :: UP(NC),VP(NC),WP(NC)		! RECONSTRUCTED VELOCITY
	REAL, INTENT(OUT) :: DIV(NC)

	INTEGER :: I
	REAL :: Q1,Q2,R1,R2,S1,S2
	REAL :: DUDX,DVDY,DWDZ
	
	! PARABOLIC EDGE RECONSTRUCTION METHOD (PERM) IN 3D, based on:
	!
	! R. McDermott and S. B. Pope. The parabolic edge reconstuction method (PERM) for Lagrangian
	! particle advection. J. Comp. Phys., 227:5447-5491, 2008.

	DO I = 1,NC 

		Q1 = Q(I)-0.5_EB
		Q2 = 0.5_EB*(Q1**2-0.25_EB)

		UP(I) = (1-S(I))*( (1-R(I))*( UBAR(1) + Q1*D1(1) + Q2*D2(1) )		&
		      +                R(I)*( UBAR(2) + Q1*D1(2) + Q2*D2(2) ) )		&
		      +     S(I)*( (1-R(I))*( UBAR(3) + Q1*D1(3) + Q2*D2(3) )		&
		      +                R(I)*( UBAR(4) + Q1*D1(4) + Q2*D2(4) ) )

		R1 = R(I)-0.5_EB
		R2 = 0.5_EB*(R1**2-0.25_EB)

		VP(I) = (1-S(I))*( (1-Q(I))*( UBAR(5) + R1*D1(5) + R2*D2(5) )		&
		      +                Q(I)*( UBAR(6) + R1*D1(6) + R2*D2(6) ) )		&
		      +     S(I)*( (1-Q(I))*( UBAR(7) + R1*D1(7) + R2*D2(7) )		&
		      +                Q(I)*( UBAR(8) + R1*D1(8) + R2*D2(8) ) )

		S1 = S(I)-0.5_EB
		S2 = 0.5_EB*(S1**2-0.25_EB)

		WP(I) = (1-Q(I))*( (1-R(I))*( UBAR(9)  + S1*D1(9)  + S2*D2(9) )		&
		      +                R(I)*( UBAR(10) + S1*D1(10) + S2*D2(10) ) )	&
		      +     Q(I)*( (1-R(I))*( UBAR(11) + S1*D1(11) + S2*D2(11) )		&
		      +                R(I)*( UBAR(12) + S1*D1(12) + S2*D2(12) ) )

	END DO


	! CHECK DIVERGENCE
	DO I = 1,NC 

		Q1 = Q(I)-0.5_EB
		
		DUDX = (1-S(I))*( (1-R(I))*( D1(1) + Q1*D2(1) )		&
		     +                R(I)*( D1(2) + Q1*D2(2) ) )	&
		     +     S(I)*( (1-R(I))*( D1(3) + Q1*D2(3) )		&
		     +                R(I)*( D1(4) + Q1*D2(4) ) )

		R1 = R(I)-0.5_EB
		
		DVDY = (1-S(I))*( (1-Q(I))*( D1(5) + R1*D2(5) )		&
		     +                Q(I)*( D1(6) + R1*D2(6) ) )	&
		     +     S(I)*( (1-Q(I))*( D1(7) + R1*D2(7) )		&
		     +                Q(I)*( D1(8) + R1*D2(8) ) )

		S1 = S(I)-0.5_EB
		
		DWDZ = (1-Q(I))*( (1-R(I))*( D1(9)  + S1*D2(9) )		&
		     +                R(I)*( D1(10) + S1*D2(10) ) )		&
		     +     Q(I)*( (1-R(I))*( D1(11) + S1*D2(11) )		&
		     +                R(I)*( D1(12) + S1*D2(12) ) )

		DIV(I) = DUDX + DVDY + DWDZ

	END DO

END SUBROUTINE PERM3D


SUBROUTINE EDGE_PARAMETERS(UBAR,D1,D2,				&
						   ULES,VLES,WLES,			&
						   UHAT,VHAT,WHAT,			&
						   DHAT,VEC_THETA)
	IMPLICIT NONE

	! DIMENSION VARIABLES

	REAL, INTENT(OUT) :: UBAR(12),D1(12),D2(12)
	REAL, INTENT(IN) :: ULES(0:1),VLES(0:1),WLES(0:1)
	REAL, INTENT(IN) :: UHAT(0:1,0:1,0:1)
	REAL, INTENT(IN) :: VHAT(0:1,0:1,0:1)
	REAL, INTENT(IN) :: WHAT(0:1,0:1,0:1)
	REAL, INTENT(IN) :: DHAT(12),VEC_THETA(8)

	! LOCAL
	REAL :: UU(0:1,0:1,0:1),VV(0:1,0:1,0:1),WW(0:1,0:1,0:1)
	REAL :: VEC_A(8),VEC_B(8)
	REAL :: DU(0:1),DV(0:1),DW(0:1)
	REAL :: DCOR(12)

	!===========================================================

	! COMPUTE FINAL VERTEX VALUES

	! U COMPONENT WEST FACE
	DU(0) = ULES(0) - 0.25_EB*( UHAT(0,0,0) + UHAT(0,0,1) + UHAT(0,1,0) + UHAT(0,1,1) )
	UU(0,0,0) = UHAT(0,0,0) + DU(0)
	UU(0,0,1) = UHAT(0,0,1) + DU(0)
	UU(0,1,0) = UHAT(0,1,0) + DU(0)
	UU(0,1,1) = UHAT(0,1,1) + DU(0)

	! U COMPONENT EAST FACE
	DU(1) = ULES(1) - 0.25_EB*( UHAT(1,0,0) + UHAT(1,0,1) + UHAT(1,1,0) + UHAT(1,1,1) )
	UU(1,0,0) = UHAT(1,0,0) + DU(1)
	UU(1,0,1) = UHAT(1,0,1) + DU(1)
	UU(1,1,0) = UHAT(1,1,0) + DU(1)
	UU(1,1,1) = UHAT(1,1,1) + DU(1)

	! V COMPONENT SOUTH FACE
	DV(0) = VLES(0) - 0.25_EB*( VHAT(0,0,0) + VHAT(0,0,1) + VHAT(1,0,0) + VHAT(1,0,1) )
	VV(0,0,0) = VHAT(0,0,0) + DV(0)
	VV(0,0,1) = VHAT(0,0,1) + DV(0)
	VV(1,0,0) = VHAT(1,0,0) + DV(0)
	VV(1,0,1) = VHAT(1,0,1) + DV(0)

	! V COMPONENT NORTH FACE
	DV(1) = VLES(1) - 0.25_EB*( VHAT(0,1,0) + VHAT(0,1,1) + VHAT(1,1,0) + VHAT(1,1,1) )
	VV(0,1,0) = VHAT(0,1,0) + DV(1)
	VV(0,1,1) = VHAT(0,1,1) + DV(1)
	VV(1,1,0) = VHAT(1,1,0) + DV(1)
	VV(1,1,1) = VHAT(1,1,1) + DV(1)

	! W COMPONENT BOTTOM FACE
	DW(0) = WLES(0) - 0.25_EB*( WHAT(0,0,0) + WHAT(0,1,0) + WHAT(1,0,0) + WHAT(1,1,0) )
	WW(0,0,0) = WHAT(0,0,0) + DW(0)
	WW(0,1,0) = WHAT(0,1,0) + DW(0)
	WW(1,0,0) = WHAT(1,0,0) + DW(0)
	WW(1,1,0) = WHAT(1,1,0) + DW(0)

	! W COMPONENT TOP FACE
	DW(1) = WLES(1) - 0.25_EB*( WHAT(0,0,1) + WHAT(0,1,1) + WHAT(1,0,1) + WHAT(1,1,1) )
	WW(0,0,1) = WHAT(0,0,1) + DW(1)
	WW(0,1,1) = WHAT(0,1,1) + DW(1)
	WW(1,0,1) = WHAT(1,0,1) + DW(1)
	WW(1,1,1) = WHAT(1,1,1) + DW(1)
	
	! PARABOLIC EDGE VELOCITIES
	UBAR(1) = 0.5_EB*( UU(1,0,0) + UU(0,0,0) )
	UBAR(2) = 0.5_EB*( UU(1,1,0) + UU(0,1,0) )
	UBAR(3) = 0.5_EB*( UU(1,0,1) + UU(0,0,1) )
	UBAR(4) = 0.5_EB*( UU(1,1,1) + UU(0,1,1) )

	UBAR(5) = 0.5_EB*( VV(0,1,0) + VV(0,0,0) )
	UBAR(6) = 0.5_EB*( VV(1,1,0) + VV(1,0,0) )
	UBAR(7) = 0.5_EB*( VV(0,1,1) + VV(0,0,1) )
	UBAR(8) = 0.5_EB*( VV(1,1,1) + VV(1,0,1) )

	UBAR(9)  = 0.5_EB*( WW(0,0,1) + WW(0,0,0) )
	UBAR(10) = 0.5_EB*( WW(0,1,1) + WW(0,1,0) )
	UBAR(11) = 0.5_EB*( WW(1,0,1) + WW(1,0,0) )
	UBAR(12) = 0.5_EB*( WW(1,1,1) + WW(1,1,0) )

	! FIRST-ORDER SLOPES
	D1(1) = UU(1,0,0) - UU(0,0,0)
	D1(2) = UU(1,1,0) - UU(0,1,0)
	D1(3) = UU(1,0,1) - UU(0,0,1)
	D1(4) = UU(1,1,1) - UU(0,1,1)

	D1(5) = VV(0,1,0) - VV(0,0,0)
	D1(6) = VV(1,1,0) - VV(1,0,0)
	D1(7) = VV(0,1,1) - VV(0,0,1)
	D1(8) = VV(1,1,1) - VV(1,0,1)

	D1(9)  = WW(0,0,1) - WW(0,0,0)
	D1(10) = WW(0,1,1) - WW(0,1,0)
	D1(11) = WW(1,0,1) - WW(1,0,0)
	D1(12) = WW(1,1,1) - WW(1,1,0)
	
	! A*D1
	VEC_A(1)  = D1(1) + D1(5) + D1(9)
	VEC_A(2)  = D1(3) + D1(7) + D1(9)
	VEC_A(3)  = D1(2) + D1(5) + D1(10)
	VEC_A(4)  = D1(4) + D1(7) + D1(10)
	VEC_A(5)  = D1(1) + D1(6) + D1(11)
	VEC_A(6)  = D1(3) + D1(8) + D1(11)
	VEC_A(7)  = D1(2) + D1(6) + D1(12)
	VEC_A(8)  = D1(4) + D1(8) + D1(12)

	! B*DHAT
	VEC_B(1)  = -DHAT(1) - DHAT(5) - DHAT(9)
	VEC_B(2)  = -DHAT(3) - DHAT(7) + DHAT(9)
	VEC_B(3)  = -DHAT(2) + DHAT(5) - DHAT(10)
	VEC_B(4)  = -DHAT(4) + DHAT(7) + DHAT(10)
	VEC_B(5)  =  DHAT(1) - DHAT(6) - DHAT(11)
	VEC_B(6)  =  DHAT(3) - DHAT(8) + DHAT(11)
	VEC_B(7)  =  DHAT(2) + DHAT(6) - DHAT(12)
	VEC_B(8)  =  DHAT(4) + DHAT(8) + DHAT(12)
	
	VEC_B = 0.5_EB*VEC_B
	
	! SECOND-ORDER SLOPES
	DCOR = MATMUL(MATRIX_BPLUS,(VEC_THETA - VEC_A - VEC_B))
	D2 = DHAT + DCOR

END SUBROUTINE EDGE_PARAMETERS


SUBROUTINE DUMP_RAW_DATA(NM)

! Write raw data files

INTEGER, INTENT(IN) :: NM
INTEGER :: LU_RAW_DATA,I,J,K
CHARACTER(80) :: FN_RAW_DATA

CALL POINT_TO_MESH(NM)

LU_RAW_DATA  = GET_FILE_NUMBER()
FN_RAW_DATA  = TRIM(CHID)//'_rawu.csv'

OPEN (UNIT=LU_RAW_DATA,FILE=FN_RAW_DATA,FORM='FORMATTED',STATUS='UNKNOWN',POSITION='REWIND')
WRITE (LU_RAW_DATA,*) '      X (m),       Y (m),       Z (m),       U (m/s)'
DO K=1,KBAR
   DO J=1,JBAR
      DO I=0,IBAR
         WRITE (LU_RAW_DATA,'(3(F12.6,A),1F14.8)') X(I),',',YC(J),',',ZC(K),',',U(I,J,K)
      ENDDO
   ENDDO
ENDDO
CLOSE (LU_RAW_DATA)

LU_RAW_DATA  = GET_FILE_NUMBER()
FN_RAW_DATA  = TRIM(CHID)//'_rawv.csv'

OPEN (UNIT=LU_RAW_DATA,FILE=FN_RAW_DATA,FORM='FORMATTED',STATUS='UNKNOWN',POSITION='REWIND')
WRITE (LU_RAW_DATA,*) '      X (m),       Y (m),       Z (m),       V (m/s)'
DO K=1,KBAR
   DO J=0,JBAR
      DO I=1,IBAR
         WRITE (LU_RAW_DATA,'(3(F12.6,A),1F14.8)') XC(I),',',Y(J),',',ZC(K),',',V(I,J,K)
      ENDDO
   ENDDO
ENDDO
CLOSE (LU_RAW_DATA)

LU_RAW_DATA  = GET_FILE_NUMBER()
FN_RAW_DATA  = TRIM(CHID)//'_raww.csv'

OPEN (UNIT=LU_RAW_DATA,FILE=FN_RAW_DATA,FORM='FORMATTED',STATUS='UNKNOWN',POSITION='REWIND')
WRITE (LU_RAW_DATA,*) '      X (m),       Y (m),       Z (m),       W (m/s)'
DO K=0,KBAR
   DO J=1,JBAR
      DO I=1,IBAR
         WRITE (LU_RAW_DATA,'(3(F12.6,A),1F14.8)') XC(I),',',YC(J),',',Z(K),',',W(I,J,K)
      ENDDO
   ENDDO
ENDDO
CLOSE (LU_RAW_DATA)

END SUBROUTINE DUMP_RAW_DATA


SUBROUTINE GET_REV_turb(MODULE_REV,MODULE_DATE)
INTEGER,INTENT(INOUT) :: MODULE_REV
CHARACTER(255),INTENT(INOUT) :: MODULE_DATE

WRITE(MODULE_DATE,'(A)') turbrev(INDEX(turbrev,':')+1:LEN_TRIM(turbrev)-2)
READ (MODULE_DATE,'(I5)') MODULE_REV
WRITE(MODULE_DATE,'(A)') turbdate

END SUBROUTINE GET_REV_turb


END MODULE TURBULENCE


! This module is an experimental implementation of my embedded mesh method (EMB),
! a prelude to adaptive mesh refinement.

MODULE EMBEDDED_MESH_METHOD

USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_VARIABLES
USE MESH_POINTERS

IMPLICIT NONE

PRIVATE
PUBLIC SCALARF_EMB,VELOCITY_EMB,RESTRICT_MASS_EMB,RESTRICT_DIV_EMB,PROJECT_VELOCITY, &
       SORT_MESH_LEVEL,MATCH_VELOCITY_EMB,SCALAR_GHOST_EMB
 
CONTAINS


SUBROUTINE SCALARF_EMB(NM1,NM2,IERROR)
IMPLICIT NONE

INTEGER, INTENT(IN) :: NM1,NM2

TYPE(MESH_TYPE), POINTER :: M1,M2
INTEGER :: N,I,J,K,I_LO,I_HI,J_LO,J_HI,K_LO,K_HI,II_0,JJ_0,KK_0,II,JJ,KK, &
           NRX,NRY,NRZ,N2X,N2Y,N2Z,II_LO,JJ_LO,KK_LO,INDEX_LIST(12),IERROR
REAL(EB) :: VOLUME_LIST(3)
REAL(EB), POINTER, DIMENSION(:,:,:,:) :: FX1,FY1,FZ1,FX2,FY2,FZ2

!   Comments:
!
!   Assumes uniform grid in each direction and that M2 lies within M1.
!
!   -------------------------------
!   |         |         |         |
!   |         |         |         |<---MESHES(M1)
!   |         |         |         |
!   |         |         |         |
!   -------------------------------
!   |         |-|-|-|-|-|         |
!   |         |-|-|-|-|-|<-------------MESHES(M2)
!   |         |-|-|-|-|-|         |
!   |         |-|-|-|-|-|         |
!   -------------------------------
!   |         |         |         |
!   |         |         |         |
!   |         |         |         |
!   |         |         |         |
!   -------------------------------

CALL LOCATE_MESH(INDEX_LIST,VOLUME_LIST,NM1,NM2,IERROR)
SELECT CASE (IERROR)
   CASE(0)
      I_LO = INDEX_LIST(1)
      I_HI = INDEX_LIST(2)
      J_LO = INDEX_LIST(3)
      J_HI = INDEX_LIST(4)
      K_LO = INDEX_LIST(5)
      K_HI = INDEX_LIST(6)
      II_LO = INDEX_LIST(7)
      JJ_LO = INDEX_LIST(8)
      KK_LO = INDEX_LIST(9)
      NRX = INDEX_LIST(10)
      NRY = INDEX_LIST(11)
      NRZ = INDEX_LIST(12)
   CASE(1)
      RETURN
END SELECT

M1=>MESHES(NM1) ! coarse mesh
M2=>MESHES(NM2) ! fine mesh

N2X = NRY*NRZ
N2Y = NRX*NRZ
N2Z = NRX*NRY

! Fluxes

FX1=>M1%SCALAR_SAVE1
FY1=>M1%SCALAR_SAVE2
FZ1=>M1%SCALAR_SAVE3

FX2=>M2%SCALAR_SAVE1
FY2=>M2%SCALAR_SAVE2
FZ2=>M2%SCALAR_SAVE3

! Restrict fine mesh to coarse mesh for embedded cells

SPECIES_LOOP: DO N=0,N_SPECIES

   ! x-direction fluxes

   DO K = K_LO,K_HI
      KK_0 = KK_LO + (K-K_LO)*NRZ
      DO J = J_LO,J_HI
         JJ_0 = JJ_LO + (J-J_LO)*NRY
         DO I = I_LO-1,I_HI !! note: this includes fine mesh boundary
            II_0 = II_LO + (I-I_LO+1)*NRX !!
                  
            FX1(I,J,K,N) = 0._EB
            DO KK = KK_0+1,KK_0+NRZ
               DO JJ = JJ_0+1,JJ_0+NRY
                  FX1(I,J,K,N) = FX1(I,J,K,N) + FX2(II_0,JJ,KK,N)
               ENDDO
            ENDDO
            FX1(I,J,K,N) = FX1(I,J,K,N)/N2X
         
         ENDDO
      ENDDO
   ENDDO
   
   ! y-direction fluxes

   DO K = K_LO,K_HI
      KK_0 = KK_LO + (K-K_LO)*NRZ
      DO J = J_LO-1,J_HI !!
         JJ_0 = JJ_LO + (J-J_LO+1)*NRY !!
         DO I = I_LO,I_HI 
            II_0 = II_LO + (I-I_LO)*NRX
                  
            FY1(I,J,K,N) = 0._EB
            DO KK = KK_0+1,KK_0+NRZ
               DO II = II_0+1,II_0+NRX
                  FY1(I,J,K,N) = FY1(I,J,K,N) + FY2(II,JJ_0,KK,N)
               ENDDO
            ENDDO
            FY1(I,J,K,N) = FY1(I,J,K,N)/N2Y
         
         ENDDO
      ENDDO
   ENDDO
   
   ! z-direction fluxes

   DO K = K_LO-1,K_HI !!
      KK_0 = KK_LO + (K-K_LO+1)*NRZ !!
      DO J = J_LO,J_HI
         JJ_0 = JJ_LO + (J-J_LO)*NRY
         DO I = I_LO,I_HI 
            II_0 = II_LO + (I-I_LO)*NRX
                  
            FZ1(I,J,K,N) = 0._EB
            DO JJ = JJ_0+1,JJ_0+NRY
               DO II = II_0+1,II_0+NRX
                  FZ1(I,J,K,N) = FZ1(I,J,K,N) + FZ2(II,JJ,KK_0,N)
               ENDDO
            ENDDO
            FZ1(I,J,K,N) = FZ1(I,J,K,N)/N2Z
         
         ENDDO
      ENDDO
   ENDDO
   
ENDDO SPECIES_LOOP

END SUBROUTINE SCALARF_EMB


SUBROUTINE VELOCITY_EMB(NM1,NM2,IERROR)
IMPLICIT NONE

INTEGER, INTENT(IN) :: NM1,NM2

TYPE(MESH_TYPE), POINTER :: M1,M2
INTEGER :: I,J,K,I_LO,I_HI,J_LO,J_HI,K_LO,K_HI,II_0,JJ_0,KK_0,II,JJ,KK, &
           NRX,NRY,NRZ,N2X,N2Y,N2Z,II_LO,JJ_LO,KK_LO,INDEX_LIST(12),IERROR
REAL(EB) :: VOLUME_LIST(3)
REAL(EB), POINTER, DIMENSION(:,:,:) :: UU1,VV1,WW1,UU2,VV2,WW2

CALL LOCATE_MESH(INDEX_LIST,VOLUME_LIST,NM1,NM2,IERROR)
SELECT CASE (IERROR)
   CASE(0)
      I_LO = INDEX_LIST(1)
      I_HI = INDEX_LIST(2)
      J_LO = INDEX_LIST(3)
      J_HI = INDEX_LIST(4)
      K_LO = INDEX_LIST(5)
      K_HI = INDEX_LIST(6)
      II_LO = INDEX_LIST(7)
      JJ_LO = INDEX_LIST(8)
      KK_LO = INDEX_LIST(9)
      NRX = INDEX_LIST(10)
      NRY = INDEX_LIST(11)
      NRZ = INDEX_LIST(12)
   CASE(1)
      RETURN
END SELECT

M1=>MESHES(NM1) ! coarse mesh
M2=>MESHES(NM2) ! fine mesh

N2X = NRY*NRZ
N2Y = NRX*NRZ
N2Z = NRX*NRY

IF (PREDICTOR) THEN
   UU1=>M1%U
   VV1=>M1%V
   WW1=>M1%W
   UU2=>M2%U
   VV2=>M2%V
   WW2=>M2%W
ELSEIF (CORRECTOR) THEN
   UU1=>M1%US
   VV1=>M1%VS
   WW1=>M1%WS
   UU2=>M2%US
   VV2=>M2%VS
   WW2=>M2%WS
ENDIF

! Restrict fine mesh to coarse mesh for embedded cells

! U-VELOCITY

DO K = K_LO,K_HI
   KK_0 = KK_LO + (K-K_LO)*NRZ
   DO J = J_LO,J_HI
      JJ_0 = JJ_LO + (J-J_LO)*NRY
      DO I = I_LO,I_HI-1 ! excludes boundary values
         II_0 = II_LO + (I-I_LO+1)*NRX
                  
         UU1(I,J,K) = 0._EB
         DO KK = KK_0+1,KK_0+NRZ
            DO JJ = JJ_0+1,JJ_0+NRY
               UU1(I,J,K) = UU1(I,J,K) + UU2(II_0,JJ,KK)
            ENDDO
         ENDDO
         UU1(I,J,K) = UU1(I,J,K)/N2X
         
      ENDDO
   ENDDO
ENDDO
   
! V-VELOCITY

DO K = K_LO,K_HI
   KK_0 = KK_LO + (K-K_LO)*NRZ
   DO J = J_LO,J_HI-1 ! excludes boundary values
      JJ_0 = JJ_LO + (J-J_LO+1)*NRY
      DO I = I_LO,I_HI 
         II_0 = II_LO + (I-I_LO)*NRX
                  
         VV1(I,J,K) = 0._EB
         DO KK = KK_0+1,KK_0+NRZ
            DO II = II_0+1,II_0+NRX
               VV1(I,J,K) = VV1(I,J,K) + VV2(II,JJ_0,KK)
            ENDDO
         ENDDO
         VV1(I,J,K) = VV1(I,J,K)/N2Y
         
      ENDDO
   ENDDO
ENDDO
   
! W-VELOCITY

DO K = K_LO,K_HI-1 ! excludes boundary values
   KK_0 = KK_LO + (K-K_LO+1)*NRZ
   DO J = J_LO,J_HI
      JJ_0 = JJ_LO + (J-J_LO)*NRY
      DO I = I_LO,I_HI 
         II_0 = II_LO + (I-I_LO)*NRX
                  
         WW1(I,J,K) = 0._EB
         DO JJ = JJ_0+1,JJ_0+NRY
            DO II = II_0+1,II_0+NRX
               WW1(I,J,K) = WW1(I,J,K) + WW2(II,JJ,KK_0)
            ENDDO
         ENDDO
         WW1(I,J,K) = WW1(I,J,K)/N2Z
         
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE VELOCITY_EMB


SUBROUTINE RESTRICT_MASS_EMB(NM1,NM2,IERROR)
IMPLICIT NONE

INTEGER, INTENT(IN) :: NM1,NM2

TYPE(MESH_TYPE), POINTER :: M1,M2
INTEGER :: N,I,J,K,I_LO,I_HI,J_LO,J_HI,K_LO,K_HI,II_0,JJ_0,KK_0,II,JJ,KK, &
           NRX,NRY,NRZ, II_LO,JJ_LO,KK_LO,INDEX_LIST(12),IERROR
REAL(EB) :: DV1,DV2,DVRAT,VOLUME_LIST(3)
REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO1,RHO2
REAL(EB), POINTER, DIMENSION(:,:,:,:) :: YY1,YY2

CALL LOCATE_MESH(INDEX_LIST,VOLUME_LIST,NM1,NM2,IERROR)
SELECT CASE (IERROR)
   CASE(0)
      I_LO = INDEX_LIST(1)
      I_HI = INDEX_LIST(2)
      J_LO = INDEX_LIST(3)
      J_HI = INDEX_LIST(4)
      K_LO = INDEX_LIST(5)
      K_HI = INDEX_LIST(6)
      II_LO = INDEX_LIST(7)
      JJ_LO = INDEX_LIST(8)
      KK_LO = INDEX_LIST(9)
      NRX = INDEX_LIST(10)
      NRY = INDEX_LIST(11)
      NRZ = INDEX_LIST(12)
      DV1 = VOLUME_LIST(1)
      DV2 = VOLUME_LIST(2)
      DVRAT = VOLUME_LIST(3)
   CASE(1)
      RETURN
END SELECT

M1=>MESHES(NM1) ! coarse mesh
M2=>MESHES(NM2) ! fine mesh

IF (PREDICTOR) THEN
   RHO1 => M1%RHOS
   RHO2 => M2%RHOS
   IF (N_SPECIES>0) YY1  => M1%YYS
   IF (N_SPECIES>0) YY2  => M2%YYS
ELSEIF (CORRECTOR) THEN
   RHO1 => M1%RHO
   RHO2 => M2%RHO
   IF (N_SPECIES>0) YY1  => M1%YY
   IF (N_SPECIES>0) YY2  => M2%YY
ENDIF

DO K = K_LO,K_HI
   KK_0 = KK_LO + (K-K_LO)*NRZ
   DO J = J_LO,J_HI
      JJ_0 = JJ_LO + (J-J_LO)*NRY
      DO I = I_LO,I_HI
         II_0 = II_LO + (I-I_LO)*NRX
            
         RHO1(I,J,K) = 0._EB
         
         DO KK = KK_0+1,KK_0+NRZ
            DO JJ = JJ_0+1,JJ_0+NRY
               DO II = II_0+1,II_0+NRX
                 
                  RHO1(I,J,K) = RHO1(I,J,K) + RHO2(II,JJ,KK)*DVRAT
                     
               ENDDO
            ENDDO
         ENDDO
      
      ENDDO
   ENDDO
ENDDO

IF (N_SPECIES>0) THEN

   SPECIES_LOOP: DO N=1,N_SPECIES
   
      DO K = K_LO,K_HI
         KK_0 = KK_LO + (K-K_LO)*NRZ
         DO J = J_LO,J_HI
            JJ_0 = JJ_LO + (J-J_LO)*NRY
            DO I = I_LO,I_HI
               II_0 = II_LO + (I-I_LO)*NRX
            
               YY1(I,J,K,N) = 0._EB
         
               DO KK = KK_0+1,KK_0+NRZ
                  DO JJ = JJ_0+1,JJ_0+NRY
                     DO II = II_0+1,II_0+NRX
                 
                        YY1(I,J,K,N) = YY1(I,J,K,N) + RHO2(II,JJ,KK)*YY2(II,JJ,KK,N)*DV2
                     
                     ENDDO
                  ENDDO
               ENDDO
               
               YY1(I,J,K,N) = YY1(I,J,K,N)/(RHO1(I,J,K)*DV1)
      
            ENDDO
         ENDDO
      ENDDO

   ENDDO SPECIES_LOOP

ENDIF

END SUBROUTINE RESTRICT_MASS_EMB


SUBROUTINE RESTRICT_DIV_EMB(NM1,NM2,IERROR)
IMPLICIT NONE

INTEGER, INTENT(IN) :: NM1,NM2

TYPE(MESH_TYPE), POINTER :: M1,M2
INTEGER :: I,J,K,I_LO,I_HI,J_LO,J_HI,K_LO,K_HI,II_0,JJ_0,KK_0,II,JJ,KK, &
           NRX,NRY,NRZ,II_LO,JJ_LO,KK_LO,INDEX_LIST(12),IERROR
REAL(EB) :: DVRAT,VOLUME_LIST(3)
REAL(EB), POINTER, DIMENSION(:,:,:) :: DP1,DP2

CALL LOCATE_MESH(INDEX_LIST,VOLUME_LIST,NM1,NM2,IERROR)
SELECT CASE (IERROR)
   CASE(0)
      I_LO = INDEX_LIST(1)
      I_HI = INDEX_LIST(2)
      J_LO = INDEX_LIST(3)
      J_HI = INDEX_LIST(4)
      K_LO = INDEX_LIST(5)
      K_HI = INDEX_LIST(6)
      II_LO = INDEX_LIST(7)
      JJ_LO = INDEX_LIST(8)
      KK_LO = INDEX_LIST(9)
      NRX = INDEX_LIST(10)
      NRY = INDEX_LIST(11)
      NRZ = INDEX_LIST(12)
      DVRAT = VOLUME_LIST(3)
   CASE(1)
      RETURN
END SELECT

M1=>MESHES(NM1) ! coarse mesh
M2=>MESHES(NM2) ! fine mesh

IF (PREDICTOR) THEN
   DP1 => M1%DS
   DP2 => M2%DS
ELSEIF (CORRECTOR) THEN
   DP1 => M1%DDDT
   DP2 => M2%DDDT
ENDIF

! Restrict divergence
   
DO K = K_LO,K_HI
   KK_0 = KK_LO + (K-K_LO)*NRZ
   DO J = J_LO,J_HI
      JJ_0 = JJ_LO + (J-J_LO)*NRY
      DO I = I_LO,I_HI
         II_0 = II_LO + (I-I_LO)*NRX
            
         DP1(I,J,K) = 0._EB
         
         DO KK = KK_0+1,KK_0+NRZ
            DO JJ = JJ_0+1,JJ_0+NRY
               DO II = II_0+1,II_0+NRX
                 
                  DP1(I,J,K) = DP1(I,J,K) + DP2(II,JJ,KK)*DVRAT
                     
               ENDDO
            ENDDO
         ENDDO
      
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE RESTRICT_DIV_EMB


SUBROUTINE PROJECT_VELOCITY(NM)
USE POIS, ONLY: H3CZSS,H2CZSS
IMPLICIT NONE

INTEGER, INTENT(IN) :: NM
INTEGER :: I,J,K
REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,DP,PP,PRHS_SAVE
REAL(EB) :: DIV,LHSS,RES,POIS_ERR

CALL POINT_TO_MESH(NM)

IF (PREDICTOR) THEN
   ! note: PROJECT_VELOCITY is called AFTER the predictor update of velocity
   UU=>US
   VV=>VS
   WW=>WS
   DP=>D
ELSEIF (CORRECTOR) THEN
   UU=>U
   VV=>V
   WW=>W
   DP=>DS
ENDIF
PP=>WORK1
PRHS_SAVE=>WORK2

! build source

DO K=1,KBAR
   DO J=1,JBAR
      DO I=1,IBAR
         DIV = (UU(I,J,K)-UU(I-1,J,K))*RDX(I) + (VV(I,J,K)-VV(I,J-1,K))*RDY(J) + (WW(I,J,K)-WW(I,J,K-1))*RDZ(K)
         PRHS(I,J,K) = DIV-DP(I,J,K)
      ENDDO
   ENDDO
ENDDO

! solve Poisson equation

BXS = 0._EB
BXF = 0._EB
BYS = 0._EB
BYF = 0._EB
BZS = 0._EB
BZF = 0._EB

PRHS_SAVE(1:IBAR,1:JBAR,1:KBAR) = PRHS
IF (.NOT.TWO_D) CALL H3CZSS(BXS,BXF,BYS,BYF,BZS,BZF,ITRN,JTRN,PRHS,POIS_PTB,SAVE2,WORK,HX)
IF (TWO_D .AND. .NOT. CYLINDRICAL) CALL H2CZSS(BXS,BXF,BZS,BZF,ITRN,PRHS,POIS_PTB,SAVE2,WORK,HX)
PP(1:IBAR,1:JBAR,1:KBAR) = PRHS

! Apply boundary conditions to PP
 
DO K=1,KBAR
   DO J=1,JBAR
      PP(0,J,K)    = -PP(1,J,K) ! use minus if Dirichlet, plus if Neumann, see init of SAVE2
      PP(IBP1,J,K) = -PP(IBAR,J,K)
   ENDDO
ENDDO
 
DO K=1,KBAR
   DO I=1,IBAR
      PP(I,0,K)    = -PP(I,1,K)
      PP(I,JBP1,K) = -PP(I,JBAR,K)
   ENDDO
ENDDO

DO J=1,JBAR
   DO I=1,IBAR
      PP(I,J,0)    = -PP(I,J,1)
      PP(I,J,KBP1) = -PP(I,J,KBAR)
   ENDDO
ENDDO

! ************************* Check the Solution *************************
 
IF (.FALSE.) THEN

   POIS_ERR = 0._EB
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            LHSS = ((PP(I+1,J,K)-PP(I,J,K))*RDXN(I) - (PP(I,J,K)-PP(I-1,J,K))*RDXN(I-1) )*RDX(I) &
                 + ((PP(I,J+1,K)-PP(I,J,K))*RDYN(J) - (PP(I,J,K)-PP(I,J-1,K))*RDYN(J-1) )*RDY(J) &
                 + ((PP(I,J,K+1)-PP(I,J,K))*RDZN(K) - (PP(I,J,K)-PP(I,J,K-1))*RDZN(K-1) )*RDZ(K)
            RES = ABS(PRHS_SAVE(I,J,K)-LHSS)
            POIS_ERR = MAX(RES,POIS_ERR)
         ENDDO
      ENDDO
   ENDDO
   WRITE(0,*) 'POIS ERROR:',pois_ptb,pois_err

ENDIF

! correct velocities

DO K=1,KBAR
   DO J=1,JBAR
      DO I=0,IBAR
         UU(I,J,K) = UU(I,J,K) - RDXN(I)*(PP(I+1,J,K)-PP(I,J,K))
      ENDDO
   ENDDO
ENDDO

DO K=1,KBAR
   DO J=0,JBAR
      DO I=1,IBAR
         VV(I,J,K) = VV(I,J,K) - RDYN(J)*(PP(I,J+1,K)-PP(I,J,K))
      ENDDO
   ENDDO
ENDDO

DO K=0,KBAR
   DO J=1,JBAR
      DO I=1,IBAR
         WW(I,J,K) = WW(I,J,K) - RDZN(K)*(PP(I,J,K+1)-PP(I,J,K))
      ENDDO
   ENDDO
ENDDO

! check divergence

IF (.FALSE.) THEN
   POIS_ERR = 0._EB
   DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            DIV = (UU(I,J,K)-UU(I-1,J,K))*RDX(I) + (VV(I,J,K)-VV(I,J-1,K))*RDY(J) + (WW(I,J,K)-WW(I,J,K-1))*RDZ(K)
            RES = ABS(DIV-DP(I,J,K))
            POIS_ERR = MAX(RES,POIS_ERR)
         ENDDO
      ENDDO
   ENDDO
   WRITE(0,*) NM,MAXVAL(ABS(DP)),POIS_ERR

ENDIF

END SUBROUTINE PROJECT_VELOCITY


SUBROUTINE SORT_MESH_LEVEL
IMPLICIT NONE

INTEGER :: IRANK,NM,ML,MLMIN,MLMAX

MESH_LIST_EMB = 0

MLMAX = MAXVAL(MESHES(1:NMESHES)%MESH_LEVEL)
MLMIN = MINVAL(MESHES(1:NMESHES)%MESH_LEVEL)

IRANK=0

DO ML=MLMAX,MLMIN,-1
   DO NM=1,NMESHES
   
      IF (MESHES(NM)%MESH_LEVEL==ML) THEN
         IRANK=IRANK+1
         MESH_LIST_EMB(IRANK)=NM
      ENDIF

   ENDDO
ENDDO

!PRINT *,MLMIN,MLMAX
!DO IRANK=1,NMESHES
!   PRINT *,MESH_LIST_EMB(IRANK)
!ENDDO
!STOP

END SUBROUTINE SORT_MESH_LEVEL


SUBROUTINE MATCH_VELOCITY_EMB(NM1,NM2,IERROR)
IMPLICIT NONE

INTEGER, INTENT(IN) :: NM1,NM2

TYPE(MESH_TYPE), POINTER :: M1,M2
INTEGER :: I,J,K,I_LO,I_HI,J_LO,J_HI,K_LO,K_HI,II_0,JJ_0,KK_0,II,JJ,KK, &
           NRX,NRY,NRZ,II_LO,JJ_LO,KK_LO,INDEX_LIST(12),IERROR,IW,IOR
REAL(EB) :: VOLUME_LIST(3)
REAL(EB), POINTER, DIMENSION(:,:,:) :: UU1,VV1,WW1,UU2,VV2,WW2
REAL(EB), PARAMETER :: RF=0.5_EB,OMRF=0.5_EB

CALL LOCATE_MESH(INDEX_LIST,VOLUME_LIST,NM1,NM2,IERROR)
SELECT CASE (IERROR)
   CASE(0)
      I_LO = INDEX_LIST(1)
      I_HI = INDEX_LIST(2)
      J_LO = INDEX_LIST(3)
      J_HI = INDEX_LIST(4)
      K_LO = INDEX_LIST(5)
      K_HI = INDEX_LIST(6)
      II_LO = INDEX_LIST(7)
      JJ_LO = INDEX_LIST(8)
      KK_LO = INDEX_LIST(9)
      NRX = INDEX_LIST(10)
      NRY = INDEX_LIST(11)
      NRZ = INDEX_LIST(12)
   CASE(1)
      RETURN
END SELECT

M1=>MESHES(NM1) ! coarse mesh
M2=>MESHES(NM2) ! fine mesh

IF (PREDICTOR) THEN
   UU1=>M1%US
   VV1=>M1%VS
   WW1=>M1%WS
   UU2=>M2%US
   VV2=>M2%VS
   WW2=>M2%WS
ELSEIF (CORRECTOR) THEN
   UU1=>M1%U
   VV1=>M1%V
   WW1=>M1%W
   UU2=>M2%U
   VV2=>M2%V
   WW2=>M2%W
ENDIF

! Set fine mesh boundary value to corresponding coarse mesh value

! U-VELOCITY

DO K = K_LO,K_HI
   KK_0 = KK_LO + (K-K_LO)*NRZ
   DO J = J_LO,J_HI
      JJ_0 = JJ_LO + (J-J_LO)*NRY

      ! east face
      I = I_HI
      II_0 = II_LO + (I-I_LO+1)*NRX
      IF (II_0==M2%IBAR) THEN
         DO KK = KK_0+1,KK_0+NRZ
            DO JJ = JJ_0+1,JJ_0+NRY
               UU2(II_0,JJ,KK) = UU1(I,J,K)
            ENDDO
         ENDDO
      ENDIF
         
      ! west face
      I = I_LO-1
      II_0 = II_LO + (I-I_LO+1)*NRX
      IF (II_0==0) THEN
         DO KK = KK_0+1,KK_0+NRZ
            DO JJ = JJ_0+1,JJ_0+NRY
               UU2(II_0,JJ,KK) = UU1(I,J,K)
            ENDDO
         ENDDO
      ENDIF
         
   ENDDO
ENDDO
   
! V-VELOCITY

DO K = K_LO,K_HI
   KK_0 = KK_LO + (K-K_LO)*NRZ
   DO I = I_LO,I_HI
      II_0 = II_LO + (I-I_LO)*NRX

      ! north face
      J = J_HI
      JJ_0 = JJ_LO + (J-J_LO+1)*NRY
      IF (JJ_0==M2%JBAR) THEN
         DO KK = KK_0+1,KK_0+NRZ
            DO II = II_0+1,II_0+NRX
               VV2(II,JJ_0,KK) = VV1(I,J,K)
            ENDDO
         ENDDO
      ENDIF
         
      ! south face
      J = J_LO-1
      JJ_0 = JJ_LO + (J-J_LO+1)*NRY
      IF (JJ_0==0) THEN
         DO KK = KK_0+1,KK_0+NRZ
            DO II = II_0+1,II_0+NRX
               VV2(II,JJ_0,KK) = VV1(I,J,K)
            ENDDO
         ENDDO
      ENDIF
         
   ENDDO
ENDDO
   
! W-VELOCITY

DO J = J_LO,J_HI
   JJ_0 = JJ_LO + (J-J_LO)*NRY
   DO I = I_LO,I_HI
      II_0 = II_LO + (I-I_LO)*NRX

      ! top face
      K = K_HI
      KK_0 = KK_LO + (K-K_LO+1)*NRZ
      IF (KK_0==M2%KBAR) THEN
         DO JJ = JJ_0+1,JJ_0+NRY
            DO II = II_0+1,II_0+NRX
               WW2(II,JJ,KK_0) = WW1(I,J,K)
            ENDDO
         ENDDO
      ENDIF
         
      ! bottom face
      K = K_LO-1
      KK_0 = KK_LO + (K-K_LO+1)*NRZ
      IF (KK_0==0) THEN
         DO JJ = JJ_0+1,JJ_0+NRY
            DO II = II_0+1,II_0+NRX
               WW2(II,JJ,KK_0) = WW1(I,J,K)
            ENDDO
         ENDDO
      ENDIF
         
   ENDDO
ENDDO

! fine mesh boundary loop

FINE_MESH_WALL_LOOP: DO IW=1,M2%NEWC
   II  = M2%IJKW(1,IW)
   JJ  = M2%IJKW(2,IW)
   KK  = M2%IJKW(3,IW)
   IOR = M2%IJKW(4,IW)
   SELECT CASE (IOR)
      CASE(1)
         M2%UVW_SAVE(IW)=UU2(0,JJ,KK)
      CASE(-1)
         M2%UVW_SAVE(IW)=UU2(M2%IBAR,JJ,KK)
      CASE(2)
         M2%UVW_SAVE(IW)=UU2(II,0,KK)
      CASE(-2)
         M2%UVW_SAVE(IW)=UU2(II,M2%JBAR,KK)
      CASE(3)
         M2%UVW_SAVE(IW)=UU2(II,JJ,0)
      CASE(-3)
         M2%UVW_SAVE(IW)=UU2(II,JJ,M2%KBAR)
   END SELECT
ENDDO FINE_MESH_WALL_LOOP

END SUBROUTINE MATCH_VELOCITY_EMB


SUBROUTINE SCALAR_GHOST_EMB(NM1,NM2,IERROR)
IMPLICIT NONE

INTEGER, INTENT(IN) :: NM1,NM2

TYPE(MESH_TYPE), POINTER :: M1,M2
INTEGER :: N,I,J,K,I_LO,I_HI,J_LO,J_HI,K_LO,K_HI,II_0,JJ_0,KK_0,II,JJ,KK, &
           NRX,NRY,NRZ,II_LO,JJ_LO,KK_LO,INDEX_LIST(12),IERROR,IW
REAL(EB) :: VOLUME_LIST(3)
REAL(EB), POINTER, DIMENSION(:,:,:) :: RHOP1,RHOP2,TMP1,TMP2,HH1,HH2
REAL(EB), POINTER, DIMENSION(:,:,:,:) :: YYP1,YYP2

CALL LOCATE_MESH(INDEX_LIST,VOLUME_LIST,NM1,NM2,IERROR)
SELECT CASE (IERROR)
   CASE(0)
      I_LO = INDEX_LIST(1)
      I_HI = INDEX_LIST(2)
      J_LO = INDEX_LIST(3)
      J_HI = INDEX_LIST(4)
      K_LO = INDEX_LIST(5)
      K_HI = INDEX_LIST(6)
      II_LO = INDEX_LIST(7)
      JJ_LO = INDEX_LIST(8)
      KK_LO = INDEX_LIST(9)
      NRX = INDEX_LIST(10)
      NRY = INDEX_LIST(11)
      NRZ = INDEX_LIST(12)
   CASE(1)
      RETURN
END SELECT

M1=>MESHES(NM1) ! coarse mesh
M2=>MESHES(NM2) ! fine mesh

IF (PREDICTOR) THEN
   RHOP1 => M1%RHOS
   YYP1  => M1%YYS
   HH1   => M1%H
   
   RHOP2 => M2%RHOS
   YYP2  => M2%YYS
   HH2   => M2%H
ELSEIF (CORRECTOR) THEN
   RHOP1 => M1%RHO
   YYP1  => M1%YY
   HH1   => M1%HS
   
   RHOP2 => M2%RHO
   YYP2  => M2%YY
   HH2   => M2%HS
ENDIF
TMP1 => M1%TMP
TMP2 => M2%TMP


! Set fine mesh boundary value to corresponding coarse mesh value

SPECIES_LOOP: DO N=1,N_SPECIES

   DO K = K_LO,K_HI
      KK_0 = KK_LO + (K-K_LO)*NRZ
      DO J = J_LO,J_HI
         JJ_0 = JJ_LO + (J-J_LO)*NRY

         ! east face
         I = I_HI+1
         II_0 = II_LO + (I-I_LO)*NRX + 1
         IF (II_0==M2%IBP1 .AND. I_HI/=M1%IBAR) THEN
         ! if I_HI==M1%IBAR then this might be an external boundary and the ghost cell value
         ! is handled in WALL_BC; similar conditions apply below
            DO KK = KK_0+1,KK_0+NRZ
               DO JJ = JJ_0+1,JJ_0+NRY
                  RHOP2(II_0,JJ,KK) = RHOP1(I,J,K)
                  TMP2(II_0,JJ,KK) = TMP1(I,J,K)
                  HH2(II_0,JJ,KK) = HH1(I,J,K)
                  IF (N_SPECIES>0) YYP2(II_0,JJ,KK,N) = YYP1(I,J,K,N)
               ENDDO
            ENDDO
         ENDIF
         
         ! west face
         I = I_LO-1
         II_0 = II_LO + (I-I_LO+1)*NRX
         IF (II_0==0 .AND. I_LO/=1) THEN
            DO KK = KK_0+1,KK_0+NRZ
               DO JJ = JJ_0+1,JJ_0+NRY
                  RHOP2(II_0,JJ,KK) = RHOP1(I,J,K)
                  TMP2(II_0,JJ,KK) = TMP1(I,J,K)
                  HH2(II_0,JJ,KK) = HH1(I,J,K)
                  IF (N_SPECIES>0) YYP2(II_0,JJ,KK,N) = YYP1(I,J,K,N)
               ENDDO
            ENDDO
         ENDIF
         
      ENDDO
   ENDDO

   DO K = K_LO,K_HI
      KK_0 = KK_LO + (K-K_LO)*NRZ
      DO I = I_LO,I_HI
         II_0 = II_LO + (I-I_LO)*NRX

         ! north face
         J = J_HI+1
         JJ_0 = JJ_LO + (J-J_LO)*NRY + 1
         IF (JJ_0==M2%JBP1 .AND. J_HI/=M1%JBAR) THEN
            DO KK = KK_0+1,KK_0+NRZ
               DO II = II_0+1,II_0+NRX
                  RHOP2(II,JJ_0,KK) = RHOP1(I,J,K)
                  TMP2(II,JJ_0,KK) = TMP1(I,J,K)
                  HH2(II,JJ_0,KK) = HH1(I,J,K)
                  IF (N_SPECIES>0) YYP2(II,JJ_0,KK,N) = YYP1(I,J,K,N)
               ENDDO
            ENDDO
         ENDIF
         
         ! south face
         J = J_LO-1
         JJ_0 = JJ_LO + (J-J_LO+1)*NRY
         IF (JJ_0==0 .AND. J_LO/=1) THEN
            DO KK = KK_0+1,KK_0+NRZ
               DO II = II_0+1,II_0+NRX
                  RHOP2(II,JJ_0,KK) = RHOP1(I,J,K)
                  TMP2(II,JJ_0,KK) = TMP1(I,J,K)
                  HH2(II,JJ_0,KK) = HH1(I,J,K)
                  IF (N_SPECIES>0) YYP2(II,JJ_0,KK,N) = YYP1(I,J,K,N)
               ENDDO
            ENDDO
         ENDIF
         
      ENDDO
   ENDDO

   DO J = J_LO,J_HI
      JJ_0 = JJ_LO + (J-J_LO)*NRY
      DO I = I_LO,I_HI
         II_0 = II_LO + (I-I_LO)*NRX

         ! top face
         K = K_HI+1
         KK_0 = KK_LO + (K-K_LO)*NRZ + 1
         IF (KK_0==M2%KBP1  .AND. K_HI/=M1%KBAR) THEN
            DO JJ = JJ_0+1,JJ_0+NRY
               DO II = II_0+1,II_0+NRX
                  RHOP2(II,JJ,KK_0) = RHOP1(I,J,K)
                  TMP2(II,JJ,KK_0) = TMP1(I,J,K)
                  HH2(II,JJ,KK_0) = HH1(I,J,K)
                  IF (N_SPECIES>0) YYP2(II,JJ,KK_0,N) = YYP1(I,J,K,N)
               ENDDO
            ENDDO
         ENDIF
         
         ! bottom face
         K = K_LO-1
         KK_0 = KK_LO + (K-K_LO+1)*NRZ
         IF (KK_0==0 .AND. K_LO/=1) THEN
            DO JJ = JJ_0+1,JJ_0+NRY
               DO II = II_0+1,II_0+NRX
                  RHOP2(II,JJ,KK_0) = RHOP1(I,J,K)
                  TMP2(II,JJ,KK_0) = TMP1(I,J,K)
                  HH2(II,JJ,KK_0) = HH1(I,J,K)
                  IF (N_SPECIES>0) YYP2(II,JJ,KK_0,N) = YYP1(I,J,K,N)
               ENDDO
            ENDDO
         ENDIF
         
      ENDDO
   ENDDO

   WALL_LOOP: DO IW=1,M2%NWC
      IF (M2%BOUNDARY_TYPE(IW)/=INTERPOLATED_BOUNDARY) CYCLE WALL_LOOP
      II = M2%IJKW(1,IW)
      JJ = M2%IJKW(2,IW)
      KK = M2%IJKW(3,IW)
      M2%RHO_F(IW)  = RHOP2(II,JJ,KK) 
      M2%YY_F(IW,N) = YYP2(II,JJ,KK,N)
   ENDDO WALL_LOOP

ENDDO SPECIES_LOOP

END SUBROUTINE SCALAR_GHOST_EMB


SUBROUTINE LOCATE_MESH(INDEX_LIST,VOLUME_LIST,NM1,NM2,IERROR)
IMPLICIT NONE

INTEGER, INTENT(IN) :: NM1,NM2

INTEGER, INTENT(OUT) :: IERROR,INDEX_LIST(12)
REAL(EB), INTENT(OUT) :: VOLUME_LIST(3)
TYPE (MESH_TYPE), POINTER :: M1,M2

INTEGER :: I_LO,I_HI,J_LO,J_HI,K_LO,K_HI,II_LO,JJ_LO,KK_LO,NRX,NRY,NRZ
REAL(EB) :: DV1,DV2,DVRAT

IERROR=0
INDEX_LIST=0
VOLUME_LIST=0._EB

M1=>MESHES(NM1) ! coarse mesh
M2=>MESHES(NM2) ! fine mesh

! Locate fine mesh within coarse mesh

I_LO = MAX(1,       NINT((M2%XS-M1%XS)/M1%DX(1))+1 )
I_HI = MIN(M1%IBAR, NINT((M2%XF-M1%XS)/M1%DX(1))   )
IF (I_LO>M1%IBAR .OR. I_HI<1) THEN ! meshes do not overlap
   IERROR=1
   RETURN
ENDIF

J_LO = MAX(1,       NINT((M2%YS-M1%YS)/M1%DY(1))+1 )
J_HI = MIN(M1%JBAR, NINT((M2%YF-M1%YS)/M1%DY(1))   )
IF (J_LO>M1%JBAR .OR. J_HI<1) THEN ! meshes do not overlap
   IERROR=1
   RETURN
ENDIF

K_LO = MAX(1,       NINT((M2%ZS-M1%ZS)/M1%DZ(1))+1 )
K_HI = MIN(M1%KBAR, NINT((M2%ZF-M1%ZS)/M1%DZ(1))   )
IF (K_LO>M1%KBAR .OR. K_HI<1) THEN ! meshes do not overlap
   IERROR=1
   RETURN
ENDIF

! Find fine mesh off-set

II_LO = MAX(0, NINT((M1%XS-M2%XS)/M2%DX(1)) )
JJ_LO = MAX(0, NINT((M1%YS-M2%YS)/M2%DY(1)) )
KK_LO = MAX(0, NINT((M1%ZS-M2%ZS)/M2%DZ(1)) )

! Compute refinment ratio in each direction

NRX = NINT(M1%DX(1)/M2%DX(1))
NRY = NINT(M1%DY(1)/M2%DY(1))
NRZ = NINT(M1%DZ(1)/M2%DZ(1))

! Cell volumes

DV1 = M1%DX(1)*M1%DY(1)*M1%DZ(1)
DV2 = M2%DX(1)*M2%DY(1)*M2%DZ(1)
DVRAT = DV2/DV1

INDEX_LIST = (/I_LO,I_HI,J_LO,J_HI,K_LO,K_HI,II_LO,JJ_LO,KK_LO,NRX,NRY,NRZ/)
VOLUME_LIST = (/DV1,DV2,DVRAT/)

END SUBROUTINE LOCATE_MESH

END MODULE EMBEDDED_MESH_METHOD


MODULE COMPLEX_GEOMETRY

USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_VARIABLES
USE MESH_POINTERS
USE TURBULENCE

IMPLICIT NONE

PRIVATE
PUBLIC :: INIT_IBM,VELTAN2D,VELTAN3D,TRILINEAR,GETX,GETU,GETGRAD
 
CONTAINS

SUBROUTINE INIT_IBM(T,NM)
USE MEMORY_FUNCTIONS, ONLY: ChkMemErr
IMPLICIT NONE
INTEGER, INTENT(IN) :: NM
REAL(EB), INTENT(IN) :: T
INTEGER :: IZERO,I,J,K,NG
TYPE (MESH_TYPE), POINTER :: M
TYPE (GEOMETRY_TYPE), POINTER :: G
REAL(EB) :: DELTA,D2,R2,RP2,XU(3),PP(3),DP,TIME,TOL=1.E-9_EB
REAL(EB) :: X_MIN,Y_MIN,Z_MIN,X_MAX,Y_MAX,Z_MAX

TIME = T
M => MESHES(NM)

! geometry loop

GEOM_LOOP: DO NG=1,N_GEOM

   G => GEOMETRY(NG)

   IF (ICYC>1 .AND. (.NOT. G%TRANSLATE) .AND. (.NOT. G%ROTATE)) CYCLE GEOM_LOOP
   
   ! acceleration (not implemented yet)
   
   G%U = G%U0
   G%V = G%V0
   G%W = G%W0
   
   ! translation (linear for now)
   
   G%X = G%X0 + G%U*TIME
   G%Y = G%Y0 + G%V*TIME
   G%Z = G%Z0 + G%W*TIME
   
   IBM_UVWMAX = MAXVAL((/ABS(G%U),ABS(G%V),ABS(G%W)/))*RDX(1)
      
   IF (TWO_D) THEN
      !DELTA = SQRT(M%DX(1)*M%DZ(1))
      DELTA = MIN(M%DX(1),M%DZ(1))
   ELSE
      !DELTA = (M%DX(1)*M%DY(1)*M%DZ(1))**ONTH
      DELTA = MIN(M%DX(1),M%DY(1),M%DZ(1))
   ENDIF
   D2 = DELTA**2

   ! find bounding box

   SELECT_SHAPE: SELECT CASE(G%ISHAPE)
      CASE(IBOX) SELECT_SHAPE
         G%X1 = G%X1 + G%U*DT
         G%X2 = G%X2 + G%U*DT
         G%Y1 = G%Y1 + G%V*DT
         G%Y2 = G%Y2 + G%V*DT
         G%Z1 = G%Z1 + G%W*DT
         G%Z2 = G%Z2 + G%W*DT
         X_MIN = G%X1
         X_MAX = G%X2
         Y_MIN = G%Y1
         Y_MAX = G%Y2
         Z_MIN = G%Z1
         Z_MAX = G%Z2
         G%HL(1) = 0.5_EB*(X_MAX-X_MIN) + TOL
         G%HL(2) = 0.5_EB*(X_MAX-X_MIN) + TOL
         G%HL(3) = 0.5_EB*(X_MAX-X_MIN) + TOL
      CASE(ISPHERE) SELECT_SHAPE
         X_MIN = G%X-G%RADIUS
         Y_MIN = G%Y-G%RADIUS
         Z_MIN = G%Z-G%RADIUS
         X_MAX = G%X+G%RADIUS
         Y_MAX = G%Y+G%RADIUS
         Z_MAX = G%Z+G%RADIUS
         R2 = G%RADIUS**2
         IBM_UVWMAX = IBM_UVWMAX + G%RADIUS*MAXVAL((/ABS(G%OMEGA_X),ABS(G%OMEGA_Y),ABS(G%OMEGA_Z)/))*RDX(1)
      CASE(ICYLINDER) SELECT_SHAPE ! assume aligned with y axis
         G%HL(1) = ABS(G%XOR-G%X)
         G%HL(2) = ABS(G%YOR-G%Y)
         G%HL(3) = ABS(G%ZOR-G%Z)
         X_MIN = G%X-G%RADIUS
         Y_MIN = G%Y-G%HL(2)
         Z_MIN = G%Z-G%RADIUS
         X_MAX = G%X+G%RADIUS
         Y_MAX = G%Y+G%HL(2)
         Z_MAX = G%Z+G%RADIUS
         R2  = G%RADIUS**2
      CASE(IPLANE)
         X_MIN = M%XS
         Y_MIN = M%YS
         Z_MIN = M%ZS
         X_MAX = M%XF
         Y_MAX = M%YF
         Z_MAX = M%ZF
         PP   = (/G%X,G%Y,G%Z/)
         G%NN = (/G%XOR,G%YOR,G%ZOR/) - PP          ! normal vector to plane
         G%NN = G%NN/SQRT(DOT_PRODUCT(G%NN,G%NN))   ! unit normal
   END SELECT SELECT_SHAPE

   G%MIN_I(NM) = M%IBAR
   G%MIN_J(NM) = M%JBAR
   G%MIN_K(NM) = M%KBAR

   IF (X_MIN>=M%XS .AND. X_MIN<=M%XF) G%MIN_I(NM) = MAX(0,FLOOR((X_MIN-M%XS)/M%DX(1))-1)
   IF (Y_MIN>=M%YS .AND. Y_MIN<=M%YF) G%MIN_J(NM) = MAX(0,FLOOR((Y_MIN-M%YS)/M%DY(1))-1)
   IF (Z_MIN>=M%ZS .AND. Z_MIN<=M%ZF) G%MIN_K(NM) = MAX(0,FLOOR((Z_MIN-M%ZS)/M%DZ(1))-1)
   
   G%MAX_I(NM) = 0
   G%MAX_J(NM) = 0
   G%MAX_K(NM) = 0

   IF (X_MAX>=M%XS .AND. X_MAX<=M%XF) G%MAX_I(NM) = MIN(M%IBAR,CEILING((X_MAX-M%XS)/M%DX(1))+1)
   IF (Y_MAX>=M%YS .AND. Y_MAX<=M%YF) G%MAX_J(NM) = MIN(M%JBAR,CEILING((Y_MAX-M%YS)/M%DY(1))+1)
   IF (Z_MAX>=M%ZS .AND. Z_MAX<=M%ZF) G%MAX_K(NM) = MIN(M%KBAR,CEILING((Z_MAX-M%ZS)/M%DZ(1))+1)
   
   IF (TWO_D) THEN
      G%MIN_J(NM)=1
      G%MAX_J(NM)=1
   ENDIF
   
   IF ( G%MAX_I(NM)<G%MIN_I(NM) .OR. &
        G%MAX_J(NM)<G%MIN_J(NM) .OR. &
        G%MAX_K(NM)<G%MIN_K(NM) ) CYCLE GEOM_LOOP
   
   IF (G%IBM_ALLOCATED) DEALLOCATE(G%U_MASK,G%V_MASK,G%W_MASK)
   G%IBM_ALLOCATED=.FALSE.
   
   ALLOCATE(G%U_MASK(G%MIN_I(NM):G%MAX_I(NM),G%MIN_J(NM):G%MAX_J(NM),G%MIN_K(NM):G%MAX_K(NM)),STAT=IZERO)
   CALL ChkMemErr('INIT_IBM','U_MASK',IZERO)
   ALLOCATE(G%V_MASK(G%MIN_I(NM):G%MAX_I(NM),G%MIN_J(NM):G%MAX_J(NM),G%MIN_K(NM):G%MAX_K(NM)),STAT=IZERO)
   CALL ChkMemErr('INIT_IBM','V_MASK',IZERO)
   ALLOCATE(G%W_MASK(G%MIN_I(NM):G%MAX_I(NM),G%MIN_J(NM):G%MAX_J(NM),G%MIN_K(NM):G%MAX_K(NM)),STAT=IZERO)
   CALL ChkMemErr('INIT_IBM','W_MASK',IZERO)
   IF (IZERO==0) G%IBM_ALLOCATED=.TRUE.
   
   G%U_MASK = 1 ! default to gas phase
   G%V_MASK = 1
   G%W_MASK = 1
   
   ! mask cells

   DO K=G%MIN_K(NM),G%MAX_K(NM)
      DO J=G%MIN_J(NM),G%MAX_J(NM)
         DO I=G%MIN_I(NM),G%MAX_I(NM)
         
            MASK_SHAPE: SELECT CASE(G%ISHAPE)
            
               CASE(IBOX) MASK_SHAPE
                  
                  ! see if the point is inside geometry
                  IF (ABS( M%X(I)-G%X)<G%HL(1) .AND. &
                      ABS(M%YC(J)-G%Y)<G%HL(2) .AND. &
                      ABS(M%ZC(K)-G%Z)<G%HL(3)) G%U_MASK(I,J,K) = -1
                  
                  IF (ABS(M%XC(I)-G%X)<G%HL(1) .AND. &
                      ABS( M%Y(J)-G%Y)<G%HL(2) .AND. &
                      ABS(M%ZC(K)-G%Z)<G%HL(3)) G%V_MASK(I,J,K) = -1
                  
                  IF (ABS(M%XC(I)-G%X)<G%HL(1) .AND. &
                      ABS(M%YC(J)-G%Y)<G%HL(2) .AND. &
                      ABS( M%Z(K)-G%Z)<G%HL(3)) G%W_MASK(I,J,K) = -1
                  
                  ! see if the point is in surface layer
                  IF (X_MAX<M%X(I) .AND. M%X(I)<X_MAX+DELTA) G%U_MASK(I,J,K) = 0
                  IF (Y_MAX<M%Y(J) .AND. M%Y(J)<Y_MAX+DELTA) G%V_MASK(I,J,K) = 0
                  IF (Z_MAX<M%Z(K) .AND. M%Z(K)<Z_MAX+DELTA) G%W_MASK(I,J,K) = 0
                  
                  IF (X_MIN-DELTA<M%X(I) .AND. M%X(I)<X_MIN) G%U_MASK(I,J,K) = 0
                  IF (Y_MIN-DELTA<M%Y(J) .AND. M%Y(J)<Y_MIN) G%V_MASK(I,J,K) = 0
                  IF (Z_MIN-DELTA<M%Z(K) .AND. M%Z(K)<Z_MIN) G%W_MASK(I,J,K) = 0
                  
               CASE(ISPHERE) MASK_SHAPE
               
                  RP2 = (M%X(I)-G%X)**2+(M%YC(J)-G%Y)**2+(M%ZC(K)-G%Z)**2
                  IF (RP2-R2 < D2 ) G%U_MASK(I,J,K) = 0
                  IF (RP2-R2 < TOL) G%U_MASK(I,J,K) = -1
                  
                  RP2 = (M%XC(I)-G%X)**2+(M%Y(J)-G%Y)**2+(M%ZC(K)-G%Z)**2
                  IF (RP2-R2 < D2 ) G%V_MASK(I,J,K) = 0
                  IF (RP2-R2 < TOL) G%V_MASK(I,J,K) = -1
                  
                  RP2 = (M%XC(I)-G%X)**2+(M%YC(J)-G%Y)**2+(M%Z(K)-G%Z)**2
                  IF (RP2-R2 < D2 ) G%W_MASK(I,J,K) = 0
                  IF (RP2-R2 < TOL) G%W_MASK(I,J,K) = -1
                  
               CASE(ICYLINDER) MASK_SHAPE ! align with y axis for now
               
                  RP2 = (M%X(I)-G%X)**2+(M%ZC(K)-G%Z)**2
                  IF (RP2-R2 < D2 ) G%U_MASK(I,J,K) = 0
                  IF (RP2-R2 < TOL) G%U_MASK(I,J,K) = -1
                  
                  RP2 = (M%XC(I)-G%X)**2+(M%ZC(K)-G%Z)**2
                  IF (RP2-R2 < D2 ) G%V_MASK(I,J,K) = 0
                  IF (RP2-R2 < TOL) G%V_MASK(I,J,K) = -1
                  
                  RP2 = (M%XC(I)-G%X)**2+(M%Z(K)-G%Z)**2
                  IF (RP2-R2 < D2 ) G%W_MASK(I,J,K) = 0
                  IF (RP2-R2 < TOL) G%W_MASK(I,J,K) = -1
                  
               CASE(IPLANE) MASK_SHAPE
               
                  ! see Section 10.3 Schneider and Eberly
                  
                  XU = (/M%X(I),M%YC(J),M%ZC(K)/)        ! point of interest
                  IF (G%TWO_SIDED) THEN
                     DP = ABS(DOT_PRODUCT(G%NN,XU-PP))   ! distance to plane
                  ELSE
                     DP = DOT_PRODUCT(G%NN,XU-PP)        ! signed distance to plane
                  ENDIF
                  IF (DP<DELTA) G%U_MASK(I,J,K) = 0
                  IF (DP<TOL)   G%U_MASK(I,J,K) = -1
                  
                  XU = (/M%XC(I),M%Y(J),M%ZC(K)/)
                  IF (G%TWO_SIDED) THEN
                     DP = ABS(DOT_PRODUCT(G%NN,XU-PP))
                  ELSE
                     DP = DOT_PRODUCT(G%NN,XU-PP)
                  ENDIF
                  IF (DP<DELTA) G%V_MASK(I,J,K) = 0
                  IF (DP<TOL)   G%V_MASK(I,J,K) = -1
                  
                  XU = (/M%XC(I),M%YC(J),M%Z(K)/)
                  IF (G%TWO_SIDED) THEN
                     DP = ABS(DOT_PRODUCT(G%NN,XU-PP))
                  ELSE
                     DP = DOT_PRODUCT(G%NN,XU-PP)
                  ENDIF
                  IF (DP<DELTA) G%W_MASK(I,J,K) = 0
                  IF (DP<TOL)   G%W_MASK(I,J,K) = -1
                  
            END SELECT MASK_SHAPE
      
         ENDDO
      ENDDO
   ENDDO

ENDDO GEOM_LOOP

END SUBROUTINE INIT_IBM


REAL(EB) FUNCTION TRILINEAR(UU,DXI,LL)
IMPLICIT NONE

REAL(EB), INTENT(IN) :: UU(0:1,0:1,0:1),DXI(3),LL(3)

! Comments:
!
! see http://local.wasp.uwa.edu.au/~pbourke/miscellaneous/interpolation/index.html
! with appropriate scaling. LL is length of side.
!
!                       UU(1,1,1)
!        z /----------/
!        ^/          / |
!        ------------  |    Particle position
!        |          |  |
!  LL(3) |   o<-----|------- DXI = [DXI(1),DXI(2),DXI(3)]
!        |          | /        
!        |          |/      Particle property at XX = TRILINEAR
!        ------------> x
!        ^
!        |
!   X0 = [0,0,0]
!
!    UU(0,0,0)
!
!===========================================================

TRILINEAR = UU(0,0,0)*(LL(1)-DXI(1))*(LL(2)-DXI(2))*(LL(3)-DXI(3)) +	&
            UU(1,0,0)*DXI(1)*(LL(2)-DXI(2))*(LL(3)-DXI(3)) +		&
            UU(0,1,0)*(LL(1)-DXI(1))*DXI(2)*(LL(3)-DXI(3)) +		&
            UU(0,0,1)*(LL(1)-DXI(1))*(LL(2)-DXI(2))*DXI(3) +		&
            UU(1,0,1)*DXI(1)*(LL(2)-DXI(2))*DXI(3) +			&
            UU(0,1,1)*(LL(1)-DXI(1))*DXI(2)*DXI(3) +			&
            UU(1,1,0)*DXI(1)*DXI(2)*(LL(3)-DXI(3)) +			&
            UU(1,1,1)*DXI(1)*DXI(2)*DXI(3)

TRILINEAR = TRILINEAR/(LL(1)*LL(2)*LL(3))

END FUNCTION TRILINEAR


SUBROUTINE GETX(XI,XU,NG)
IMPLICIT NONE

REAL(EB), INTENT(OUT) :: XI(3)
REAL(EB), INTENT(IN) :: XU(3)
INTEGER, INTENT(IN) :: NG
TYPE(GEOMETRY_TYPE), POINTER :: G
REAL(EB) :: PP(3),RU(3),RUMAG,DR,DP

G => GEOMETRY(NG)
SELECT CASE(G%ISHAPE)
   CASE(IBOX)
      XI = XU
      IF (XU(1)<G%X1) XI(1) = XU(1) + (XU(1)-G%X1)
      IF (XU(1)>G%X2) XI(1) = XU(1) + (XU(1)-G%X2)
      IF (XU(2)<G%Y1) XI(2) = XU(2) + (XU(2)-G%Y1)
      IF (XU(2)>G%Y2) XI(2) = XU(2) + (XU(2)-G%Y2)
      IF (XU(3)<G%Z1) XI(3) = XU(3) + (XU(3)-G%Z1)
      IF (XU(3)>G%Z2) XI(3) = XU(3) + (XU(3)-G%Z2)
   CASE(ISPHERE)
      RU     = XU-(/G%X,G%Y,G%Z/)
      RUMAG  = SQRT(DOT_PRODUCT(RU,RU))
      DR     = RUMAG-G%RADIUS
      XI     = XU + DR*RU/RUMAG
   CASE(ICYLINDER)
      RU     = (/XU(1),0._EB,XU(3)/)-(/G%X,0._EB,G%Z/)
      RUMAG  = SQRT(DOT_PRODUCT(RU,RU))
      DR     = RUMAG-G%RADIUS
      XI     = XU + DR*RU/RUMAG
   CASE(IPLANE)
      PP = (/G%X,G%Y,G%Z/)           ! point in the plane
      DP = DOT_PRODUCT(G%NN,XU-PP)   ! signed distance to plane
      XI = XU + DP*G%NN
END SELECT

END SUBROUTINE GETX


SUBROUTINE GETU(U_DATA,DXI,XI,XU,INDU,I_VEL,NM)
IMPLICIT NONE

REAL(EB), INTENT(OUT) :: U_DATA(0:1,0:1,0:1),DXI(3)
REAL(EB), INTENT(IN) :: XI(3),XU(3)
INTEGER, INTENT(IN) :: INDU(3),I_VEL,NM
TYPE(MESH_TYPE), POINTER :: M
REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW
INTEGER :: II,JJ,KK

M=>MESHES(NM)
IF (PREDICTOR) THEN
   UU => M%U
   VV => M%V
   WW => M%W
ELSE
   UU => M%US
   VV => M%VS
   WW => M%WS
ENDIF

! first assume XI >= XU
II = INDU(1)
JJ = INDU(2)
KK = INDU(3)
DXI = XI-XU

SELECT CASE(I_VEL)
   CASE(1)
      IF (XI(1)<XU(1)) THEN
         II=MAX(0,II-1)
         DXI(1)=DXI(1)+M%DX(II)
      ENDIF
      IF (XI(2)<XU(2)) THEN
         JJ=MAX(0,JJ-1)
         DXI(2)=DXI(2)+M%DYN(JJ)
      ENDIF
      IF (XI(3)<XU(3)) THEN
         KK=MAX(0,KK-1)
         DXI(3)=DXI(3)+M%DZN(KK)
      ENDIF
      U_DATA(0,0,0) = UU(II,JJ,KK)
      U_DATA(1,0,0) = UU(II+1,JJ,KK)
      U_DATA(0,1,0) = UU(II,JJ+1,KK)
      U_DATA(0,0,1) = UU(II,JJ,KK+1)
      U_DATA(1,0,1) = UU(II+1,JJ,KK+1)
      U_DATA(0,1,1) = UU(II,JJ+1,KK+1)
      U_DATA(1,1,0) = UU(II+1,JJ+1,KK)
      U_DATA(1,1,1) = UU(II+1,JJ+1,KK+1)
   CASE(2)
      IF (XI(1)<XU(1)) THEN
         II=MAX(0,II-1)
         DXI(1)=DXI(1)+M%DXN(II)
      ENDIF
      IF (XI(2)<XU(2)) THEN
         JJ=MAX(0,JJ-1)
         DXI(2)=DXI(2)+M%DY(JJ)
      ENDIF
      IF (XI(3)<XU(3)) THEN
         KK=MAX(0,KK-1)
         DXI(3)=DXI(3)+M%DZN(KK)
      ENDIF
      U_DATA(0,0,0) = VV(II,JJ,KK)
      U_DATA(1,0,0) = VV(II+1,JJ,KK)
      U_DATA(0,1,0) = VV(II,JJ+1,KK)
      U_DATA(0,0,1) = VV(II,JJ,KK+1)
      U_DATA(1,0,1) = VV(II+1,JJ,KK+1)
      U_DATA(0,1,1) = VV(II,JJ+1,KK+1)
      U_DATA(1,1,0) = VV(II+1,JJ+1,KK)
      U_DATA(1,1,1) = VV(II+1,JJ+1,KK+1)
   CASE(3)
      IF (XI(1)<XU(1)) THEN
         II=MAX(0,II-1)
         DXI(1)=DXI(1)+M%DXN(II)
      ENDIF
      IF (XI(2)<XU(2)) THEN
         JJ=MAX(0,JJ-1)
         DXI(2)=DXI(2)+M%DYN(JJ)
      ENDIF
      IF (XI(3)<XU(3)) THEN
         KK=MAX(0,KK-1)
         DXI(3)=DXI(3)+M%DZ(KK)
      ENDIF
      U_DATA(0,0,0) = WW(II,JJ,KK)
      U_DATA(1,0,0) = WW(II+1,JJ,KK)
      U_DATA(0,1,0) = WW(II,JJ+1,KK)
      U_DATA(0,0,1) = WW(II,JJ,KK+1)
      U_DATA(1,0,1) = WW(II+1,JJ,KK+1)
      U_DATA(0,1,1) = WW(II,JJ+1,KK+1)
      U_DATA(1,1,0) = WW(II+1,JJ+1,KK)
      U_DATA(1,1,1) = WW(II+1,JJ+1,KK+1)
   CASE(4) ! viscosity
      IF (XI(1)<XU(1)) THEN
         II=MAX(0,II-1)
         DXI(1)=DXI(1)+M%DXN(II)
      ENDIF
      IF (XI(2)<XU(2)) THEN
         JJ=MAX(0,JJ-1)
         DXI(2)=DXI(2)+M%DYN(JJ)
      ENDIF
      IF (XI(3)<XU(3)) THEN
         KK=MAX(0,KK-1)
         DXI(3)=DXI(3)+M%DZN(KK)
      ENDIF
      U_DATA(0,0,0) = M%MU(II,JJ,KK)
      U_DATA(1,0,0) = M%MU(II+1,JJ,KK)
      U_DATA(0,1,0) = M%MU(II,JJ+1,KK)
      U_DATA(0,0,1) = M%MU(II,JJ,KK+1)
      U_DATA(1,0,1) = M%MU(II+1,JJ,KK+1)
      U_DATA(0,1,1) = M%MU(II,JJ+1,KK+1)
      U_DATA(1,1,0) = M%MU(II+1,JJ+1,KK)
      U_DATA(1,1,1) = M%MU(II+1,JJ+1,KK+1)
END SELECT

END SUBROUTINE GETU


SUBROUTINE GETGRAD(G_DATA,DXI,XI,XU,INDU,COMP_I,COMP_J,NM)
IMPLICIT NONE

REAL(EB), INTENT(OUT) :: G_DATA(0:1,0:1,0:1),DXI(3)
REAL(EB), INTENT(IN) :: XI(3),XU(3)
INTEGER, INTENT(IN) :: INDU(3),COMP_I,COMP_J,NM
TYPE(MESH_TYPE), POINTER :: M
REAL(EB), POINTER, DIMENSION(:,:,:) :: DUDX
INTEGER :: II,JJ,KK

M=>MESHES(NM)

IF (COMP_I==1 .AND. COMP_J==1) DUDX => M%WORK5
IF (COMP_I==1 .AND. COMP_J==2) DUDX => M%IBM_SAVE1
IF (COMP_I==1 .AND. COMP_J==3) DUDX => M%IBM_SAVE2
IF (COMP_I==2 .AND. COMP_J==1) DUDX => M%IBM_SAVE3
IF (COMP_I==2 .AND. COMP_J==2) DUDX => M%WORK6
IF (COMP_I==2 .AND. COMP_J==3) DUDX => M%IBM_SAVE4
IF (COMP_I==3 .AND. COMP_J==1) DUDX => M%IBM_SAVE5
IF (COMP_I==3 .AND. COMP_J==2) DUDX => M%IBM_SAVE6
IF (COMP_I==3 .AND. COMP_J==3) DUDX => M%WORK7

! first assume XI >= XU
II = INDU(1)
JJ = INDU(2)
KK = INDU(3)
DXI = XI-XU

IF (XI(1)<XU(1)) THEN
   II=MAX(0,II-1)
   DXI(1)=DXI(1)+M%DX(II)
ENDIF
IF (XI(2)<XU(2)) THEN
   JJ=MAX(0,JJ-1)
   DXI(2)=DXI(2)+M%DY(JJ)
ENDIF
IF (XI(3)<XU(3)) THEN
   KK=MAX(0,KK-1)
   DXI(3)=DXI(3)+M%DZ(KK)
ENDIF
G_DATA(0,0,0) = DUDX(II,JJ,KK)
G_DATA(1,0,0) = DUDX(II+1,JJ,KK)
G_DATA(0,1,0) = DUDX(II,JJ+1,KK)
G_DATA(0,0,1) = DUDX(II,JJ,KK+1)
G_DATA(1,0,1) = DUDX(II+1,JJ,KK+1)
G_DATA(0,1,1) = DUDX(II,JJ+1,KK+1)
G_DATA(1,1,0) = DUDX(II+1,JJ+1,KK)
G_DATA(1,1,1) = DUDX(II+1,JJ+1,KK+1)

END SUBROUTINE GETGRAD


REAL(EB) FUNCTION VELTAN2D(UU,UW,NN,DN,DIVU,GRADU,GRADP,TAU_IJ,DT,RRHO,MU,I_VEL)
IMPLICIT NONE

REAL(EB), INTENT(IN) :: UU(2),UW(2),NN(2),DN,DIVU,GRADU(2,2),GRADP(2),TAU_IJ(2,2),DT,RRHO,MU
INTEGER, INTENT(IN) :: I_VEL
REAL(EB) :: C(2,2),SS(2),SLIP_COEF,ETA,AA,BB,US0,DUMMY
REAL(EB) :: US,UN,US_WALL,UN_WALL,DPDS,DUSDS,DUSDN,TSN
INTEGER :: SUBIT

! Cartesian grid coordinate system orthonormal basis vectors
REAL(EB), DIMENSION(2), PARAMETER :: XX=(/1._EB, 0._EB/),YY=(/0._EB, 1._EB/)


! streamwise unit vector
SS = (/NN(2),-NN(1)/)

! directional cosines (see Pope, Eq. A.11)
C(1,1) = DOT_PRODUCT(XX,SS)
C(1,2) = DOT_PRODUCT(XX,NN)
C(2,1) = DOT_PRODUCT(YY,SS)
C(2,2) = DOT_PRODUCT(YY,NN)

! transform velocity (see Pope, Eq. A.17)
US = C(1,1)*UU(1) + C(2,1)*UU(2)
UN = C(1,2)*UU(1) + C(2,2)*UU(2)

! transform wall velocity
US_WALL = C(1,1)*UW(1) + C(2,1)*UW(2)
UN_WALL = C(1,2)*UW(1) + C(2,2)*UW(2)

! transform pressure gradient
DPDS = C(1,1)*GRADP(1) + C(2,1)*GRADP(2)

! transform velocity gradient tensor (Pope A.23)
DUSDS = C(1,1)*C(1,1)*GRADU(1,1) + C(1,1)*C(2,1)*GRADU(1,2) &
      + C(2,1)*C(1,1)*GRADU(2,1) + C(2,1)*C(2,1)*GRADU(2,2)
      
DUSDN = C(1,1)*C(1,2)*GRADU(1,1) + C(1,1)*C(2,2)*GRADU(1,2) &
      + C(2,1)*C(1,2)*GRADU(2,1) + C(2,1)*C(2,2)*GRADU(2,2)

! transform stress tensor
TSN = C(1,1)*C(1,2)*TAU_IJ(1,1) + C(1,1)*C(2,2)*TAU_IJ(1,2) &
    + C(2,1)*C(1,2)*TAU_IJ(2,1) + C(2,1)*C(2,2)*TAU_IJ(2,2)
    
! update boundary layer equations

! update wall-normal velocity
UN = UN_WALL + DN*(DIVU-0.5_EB*DUSDS)

! ODE solution
IF (DNS) THEN
   ETA = UN + RRHO*MU/DN
   AA  = -(0.5_EB*DUSDS + TWTH*ETA/DN)
   BB  = (TWTH*US_WALL/DN + ONSI*DUSDN)*ETA - (UN*0.5_EB*DUSDN + RRHO*( DPDS + TSN/(2._EB*DN) ))
   !AA  = -0.5_EB*(DUSDS + ETA/DN)
   !BB  = 0.5_EB*US_WALL/DN*ETA - (UN*0.5_EB*DUSDN + RRHO*( DPDS + TSN/(2._EB*DN) ))
   US = ((AA*US + BB)*EXP(AA*DT) - BB)/AA
ELSE
   US0 = US
   DO SUBIT=1,1
      CALL WERNER_WENGLE_WALL_MODEL(SLIP_COEF,DUMMY,US-US_WALL,MU*RRHO,DN,0._EB)
      !IF (SLIP_COEF< -1._EB .OR. SLIP_COEF>-1._EB) THEN
      !   PRINT *,SUBIT,'WARNING: SLIP_COEF=',SLIP_COEF
      !ENDIF
      ETA = RRHO*(1-SLIP_COEF)*MU/(2._EB*DN**2)
      AA  = -(0.5_EB*DUSDS + TWTH*UN/DN + ETA)
      BB  = ETA*US_WALL - (UN*ONTH*DUSDN + RRHO*( DPDS + TSN/(2._EB*DN) ))
      US = ((AA*US0 + BB)*EXP(AA*DT) - BB)/AA
   ENDDO
ENDIF

! transform velocity back to Cartesian component I_VEL
VELTAN2D = C(I_VEL,1)*US + C(I_VEL,2)*UN

END FUNCTION VELTAN2D


REAL(EB) FUNCTION VELTAN3D(UU,UW,NN,DN,DIVU,GRADU,GRADP,TAU_IJ,DT,RRHO,MU,I_VEL,ROUGHNESS)
USE MATH_FUNCTIONS, ONLY: CROSS_PRODUCT, NORM2
IMPLICIT NONE

REAL(EB), INTENT(IN) :: UU(3),UW(3),NN(3),DN,DIVU,GRADU(3,3),GRADP(3),TAU_IJ(3,3),DT,RRHO,MU,ROUGHNESS
INTEGER, INTENT(IN) :: I_VEL
REAL(EB) :: C(3,3),SS(3),PP(3),SLIP_COEF,ETA,AA,BB,US0,DUMMY,UU_REL(3)
REAL(EB) :: US,UP,UN,DPDS,DUSDS,DUSDN,TSN
INTEGER :: SUBIT,I,J

! Cartesian grid coordinate system orthonormal basis vectors
REAL(EB), DIMENSION(3), PARAMETER :: E1=(/1._EB,0._EB,0._EB/),E2=(/0._EB,1._EB,0._EB/),E3=(/0._EB,0._EB,1._EB/)

UU_REL = UU-UW
IF (DOT_PRODUCT(NN,UU_REL)<1.E-6_EB) THEN
   VELTAN3D = UU(I_VEL)
   RETURN
ENDIF

! find a vector PP in the tangent plane of the surface and orthogonal to UU-UW
CALL CROSS_PRODUCT(PP,NN,UU_REL)
PP = PP/NORM2(PP) ! normalize to unit vector

! define the streamwise unit vector SS
CALL CROSS_PRODUCT(SS,PP,NN)

! directional cosines (see Pope, Eq. A.11)
C(1,1) = DOT_PRODUCT(E1,SS)
C(1,2) = DOT_PRODUCT(E1,PP)
C(1,3) = DOT_PRODUCT(E1,NN)
C(2,1) = DOT_PRODUCT(E2,SS)
C(2,2) = DOT_PRODUCT(E2,PP)
C(2,3) = DOT_PRODUCT(E2,NN)
C(3,1) = DOT_PRODUCT(E3,SS)
C(3,2) = DOT_PRODUCT(E3,PP)
C(3,3) = DOT_PRODUCT(E3,NN)

! transform velocity (see Pope, Eq. A.17)
US = C(1,1)*(UU(1)-UW(1)) + C(2,1)*(UU(2)-UW(2)) + C(3,1)*(UU(3)-UW(3))
UP = C(1,2)*(UU(1)-UW(1)) + C(2,2)*(UU(2)-UW(2)) + C(3,2)*(UU(3)-UW(3))
UN = C(1,3)*(UU(1)-UW(1)) + C(2,3)*(UU(2)-UW(2)) + C(3,3)*(UU(3)-UW(3))

!! check UP, should be zero
!print *, UP

! transform pressure gradient
DPDS = C(1,1)*GRADP(1) + C(2,1)*GRADP(2) + C(3,1)*GRADP(3)

! transform tensors (Pope A.23)
DUSDS = 0._EB
DUSDN = 0._EB
TSN = 0._EB
DO I=1,3
   DO J=1,3
      DUSDS = DUSDS + C(I,1)*C(J,1)*GRADU(I,J)
      DUSDN = DUSDN + C(I,1)*C(J,3)*GRADU(I,J)
      TSN = TSN + C(I,1)*C(J,3)*TAU_IJ(I,J)
   ENDDO
ENDDO
    
! update boundary layer equations

! update wall-normal velocity
UN = DN*(DIVU-0.5_EB*DUSDS)

! ODE solution
IF (DNS) THEN
   ETA = UN + RRHO*MU/DN
   AA  = -(0.5_EB*DUSDS + TWTH*ETA/DN)
   BB  = ONSI*DUSDN*ETA - (UN*0.5_EB*DUSDN + RRHO*( DPDS + TSN/(2._EB*DN) ))
   !AA  = -0.5_EB*(DUSDS + ETA/DN)
   !BB  = 0.5_EB*US_WALL/DN*ETA - (UN*0.5_EB*DUSDN + RRHO*( DPDS + TSN/(2._EB*DN) ))
   US = ((AA*US + BB)*EXP(AA*DT) - BB)/AA
ELSE
   US0 = US
   DO SUBIT=1,1
      CALL WERNER_WENGLE_WALL_MODEL(SLIP_COEF,DUMMY,US,MU*RRHO,DN,ROUGHNESS)
      !IF (SLIP_COEF< -1._EB .OR. SLIP_COEF>-1._EB) THEN
      !   PRINT *,SUBIT,'WARNING: SLIP_COEF=',SLIP_COEF
      !ENDIF
      ETA = RRHO*(1-SLIP_COEF)*MU/(2._EB*DN**2)
      AA  = -(0.5_EB*DUSDS + TWTH*UN/DN + ETA)
      BB  = -(UN*ONTH*DUSDN + RRHO*( DPDS + TSN/(2._EB*DN) ))
      US = ((AA*US0 + BB)*EXP(AA*DT) - BB)/AA
   ENDDO
ENDIF

! transform velocity back to Cartesian component I_VEL
VELTAN3D = C(I_VEL,1)*US + C(I_VEL,3)*UN + UW(I_VEL)

END FUNCTION VELTAN3D



END MODULE COMPLEX_GEOMETRY
